Statistical Analysis of compositional data – FlareGas – 16s

16S analysis

sample.biom.silva <- paste0(devdir, "data/16s.emu.all_json_md.biom")
sample.tree.silva <- paste0(devdir, "data/16s.emu.tree.nwk")
biom.silva <- import_biom(sample.biom.silva, treefilename = sample.tree.silva)
colnames(tax_table(biom.silva)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
sampleid.silva <- sample_names(biom.silva)
labels.biom.silva <- sample_data(biom.silva)$Label

rank_names(biom.silva)

## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"

sample_variables(biom.silva)

## [1] "BarcodeSequence" "Culture_setup"   "Inoculum"        "Label"          
## [5] "PCR_amplicon"    "Primer"          "ProjectID"       "Time_point"

sample_names(biom.silva)

##  [1] "flaregas-1"  "flaregas-2"  "flaregas-3"  "flaregas-4"  "flaregas-5" 
##  [6] "flaregas-6"  "flaregas-7"  "flaregas-8"  "flaregas-9"  "flaregas-10"
## [11] "flaregas-11" "flaregas-12" "flaregas-13" "flaregas-14" "flaregas-15"
## [16] "flaregas-16" "flaregas-17" "flaregas-18" "flaregas-19" "flaregas-20"
## [21] "flaregas-21" "flaregas-22" "flaregas-23" "flaregas-24" "flaregas-25"
## [26] "flaregas-26" "flaregas-27" "flaregas-28" "flaregas-29" "flaregas-30"

knitr::kable(as.matrix(sample_data(biom.silva)))

	BarcodeSequence	Culture_setup	Inoculum	Label	PCR_amplicon	Primer	ProjectID	Time_point
flaregas-1	BC01	ocean	Helsingor	H-LC-t0	16S	27F+1492R	FAS64881	0
flaregas-2	BC02	light+clean_gas	Helsingor	H-LC-t1	16S	27F+1492R	FAS64881	7
flaregas-3	BC03	light+clean_gas	Helsingor	H-LC-t3	16S	27F+1492R	FAS64881	14
flaregas-4	BC04	ocean	Helsingor	H-LC-t0	18S	ss5+ss3	FAS64881	0
flaregas-5	BC05	light+clean_gas	Helsingor	H-LC-t1	18S	ss5+ss3	FAS64881	7
flaregas-6	BC06	light+clean_gas	Helsingor	H-LC-t3	18S	ss5+ss3	FAS64881	14
flaregas-7	BC07	ocean	Helsingor	H-LC-t0	18S	EukA+EukB	FAS64881	0
flaregas-8	BC08	light+clean_gas	Helsingor	H-LC-t1	18S	EukA+EukB	FAS64881	7
flaregas-9	BC09	light+clean_gas	Helsingor	H-LC-t3	18S	EukA+EukB	FAS64881	14
flaregas-10	BC10	ocean	Helsingor	H-LC-t0	16S+18S	27F+1492R++EukA+EukB	FAS64881	0
flaregas-11	BC01	ocean	Hornbaek-deep	Deep-DC-t0	16S+18S	27F+1492R++EukA+EukB	FAS65981	0
flaregas-12	BC02	ocean	Hornbaek-deep	Deep-DC-t0	16S+18S	27F+1492R++EukA+EukB	FAS65981	0
flaregas-13	BC03	light+clean_gas	Hornbaek-deep	Deep-LC-t1	16S+18S	27F+1492R++EukA+EukB	FAS65981	5
flaregas-14	BC04	light+clean_gas	Hornbaek-deep	Deep-LC-t2	16S+18S	27F+1492R++EukA+EukB	FAS65981	7
flaregas-15	BC05	light+clean_gas	Hornbaek-deep	Deep-LC-t3	16S+18S	27F+1492R++EukA+EukB	FAS65981	11
flaregas-16	BC06	light+clean_gas	Hornbaek-deep	Deep-LC-t3	16S+18S	27F+1492R++EukA+EukB	FAS65981	11
flaregas-17	BC07	light+used_gas	Hornbaek-deep	Deep-LU-t1	16S+18S	27F+1492R++EukA+EukB	FAS65981	5
flaregas-18	BC08	light+used_gas	Hornbaek-deep	Deep-LU-t2	16S+18S	27F+1492R++EukA+EukB	FAS65981	7
flaregas-19	BC09	light+used_gas	Hornbaek-deep	Deep-LU-t3	16S+18S	27F+1492R++EukA+EukB	FAS65981	11
flaregas-20	BC10	light+used_gas	Hornbaek-deep	Deep-LU-t3	16S+18S	27F+1492R++EukA+EukB	FAS65981	11
flaregas-21	BC11	dark+clean_gas	Hornbaek-deep	Deep-DC-t1	16S+18S	27F+1492R++EukA+EukB	FAS64830	5
flaregas-22	BC12	dark+clean_gas	Hornbaek-deep	Deep-DC-t2	16S+18S	27F+1492R++EukA+EukB	FAS64830	7
flaregas-23	BC13	dark+clean_gas	Hornbaek-deep	Deep-DC-t3	16S+18S	27F+1492R++EukA+EukB	FAS64830	11
flaregas-24	BC14	dark+clean_gas	Hornbaek-deep	Deep-DC-t3	16S+18S	27F+1492R++EukA+EukB	FAS64830	11
flaregas-25	BC15	ocean	Hornbaek-shallow	Shallow-DC-t0	16S+18S	27F+1492R++EukA+EukB	FAS64830	0
flaregas-26	BC16	ocean	Hornbaek-shallow	Shallow-DC-t0	16S+18S	27F+1492R++EukA+EukB	FAS64830	0
flaregas-27	BC17	dark+clean_gas	Hornbaek-shallow	Shallow-DC-t1	16S+18S	27F+1492R++EukA+EukB	FAS64830	5
flaregas-28	BC18	dark+clean_gas	Hornbaek-shallow	Shallow-DC-t2	16S+18S	27F+1492R++EukA+EukB	FAS64830	7
flaregas-29	BC19	dark+clean_gas	Hornbaek-shallow	Shallow-DC-t3	16S+18S	27F+1492R++EukA+EukB	FAS64830	11
flaregas-30	BC20	dark+clean_gas	Hornbaek-shallow	Shallow-DC-t3	16S+18S	27F+1492R++EukA+EukB	FAS64830	11

Preprocessing data

Subset prokaryotes (16S)

#remove Chloroplast (order level) and Mitochondria (family level)
biom.filter.silva <- subset_taxa(biom.silva, Order!="Chloroplast")
biom.filter.silva <- subset_taxa(biom.filter.silva, Family!="Mitochondria")

#get prokaryotes (16S)
biom.16s <- subset_taxa(biom.filter.silva, Kingdom=="Bacteria")

#remove 18s specific samples
biom.16s <- prune_samples(!(sample_names(biom.16s)) %in% c("flaregas-4","flaregas-5","flaregas-6","flaregas-7","flaregas-8","flaregas-9"), biom.16s)
sampleid.16s <- sample_names(biom.16s)
label.16s <- sample_data(biom.16s)$Label

meta <- rownames_to_column(as.data.frame(as.matrix(biom.16s@sam_data))) %>% dplyr::rename("Sample"="rowname")

##convert counts to integers for alpha diversity measures
biom.16s.int <- biom.16s
otu_table(biom.16s.int) <- otu_table(matrix(as.integer(otu_table(biom.16s)), ncol = nsamples(biom.16s), dimnames = list(taxa_names(biom.16s), sample_names(biom.16s))), taxa_are_rows = taxa_are_rows(biom.16s))

Define Methanotrophs

#define methanotrophs
#https://doi.org/10.1007/978-3-319-74866-5_2
methanotrophs.family <- c("Methylococcaceae",
                   "Methylohalobiaceae",
                   "Methylomonadaceae")
methanotrophs.genus <- c("Methylobacterium-Methylorubrum",
                   "Methylocapsa",
                   "Methylocella",
                   "Methylocystis",
                   "Methyloferula",
                   "Methylosinus",
                   "Methylorosula",
                   "Methylovirgula",
                   "Candidatus_Methylacidiphilum",
                   "Methyloceanibacter")

tax <- data.frame(tax_table(biom.16s))
tax.clean <- tax %>% mutate(Phylum = if_else(Genus %in% methanotrophs.genus, "Methanotroph", Phylum))
tax.clean <- tax.clean %>% mutate(Phylum = if_else(Family %in% methanotrophs.family, "Methanotroph", Phylum))
tax_table(biom.16s) <- as.matrix(tax.clean)

Prevalence filtering

Remove all taxa that only occur in a single sample

biom.16s

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1584 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 1584 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1584 tips and 1582 internal nodes ]

#tax_stats(biom.16s)
# Define prevalence of each taxa
# (in how many samples did each taxa appear at least once)
prev0 = apply(X = otu_table(biom.16s),
              MARGIN = ifelse(taxa_are_rows(biom.16s), yes = 1, no = 2),
              FUN = function(x){sum(x > 0)})
prevdf = data.frame(Prevalence = prev0, TotalAbundance = taxa_sums(biom.16s), tax_table(biom.16s))
ggplot(prevdf, aes(Prevalence)) + geom_bar()

# Define prevalence threshold as 5% of total samples
prevalenceThreshold = 0.05 * nsamples(biom.16s)
prevalenceThreshold

## [1] 1.2

# Execute prevalence filter, using `prune_taxa()` function
biom.16s = prune_taxa((prev0 > prevalenceThreshold), biom.16s)
biom.16s

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 978 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 978 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 978 tips and 976 internal nodes ]

#tax_stats(biom.16s)

#pseq <- core(pseq.comp, detection = 1/1000, prevalence = 6/30) #table(meta$Culture_setup) -> minimum no. samples of dark-light-ocean: 6
# Merge rare taxa to speed up examples
#pseq <- aggregate_rare(pseq.comp, level = "Genus", detection = 0, prevalence = 6/30)

Normalization and transformation

Relative abundance

# Transform to relative abundance. Save as new object.
biom.16s.rel <- biom.16s %>% transform_sample_counts(function(x) {x / sum(x)})
#alternatively run: biom.16s.rel <- microbiome::transform(biom.16s, "compositional")

# visualize relative abundance distribution
d.rel.melt <- melt(otu_table(biom.16s.rel), id.vars = colnames(otu_table(biom.16s.rel))) %>% dplyr::rename("Sample" = "Var2", "Relative abundance" = "value")
d.rel.melt <- inner_join(d.rel.melt, meta, by = "Sample")
pdf(paste0(wdir,"fig/rel_abundance-distr-16s.pdf"), width = 14, height = 10)
ggplot(d.rel.melt, aes(`Relative abundance`, Sample, fill=Inoculum)) + geom_boxplot(outliers = F) + coord_flip() + facet_wrap(~ Culture_setup, ncol = 2, scales = "free_x") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
dev.off()

Centered log-ratio transformation

# load the library zCompositions to perform 0 replacement
library(zCompositions)

# we are using the Count Zero Multiplicative approach
# z.warning is threshold used to identify individual rows or columns including an excess of zeros/unobserved values
d.n0 <- cmultRepl(otu_table(biom.16s), method="CZM", label=0, z.warning=0.99)
dim(d.n0)

# generate the centered log-ratio transformed data
# samples by row
d.clr <- apply(d.n0, 2, function(x) log(x) - mean(log(x)))
biom.16s.clr <- biom.16s
otu_table(biom.16s.clr) <- otu_table(d.clr, taxa_are_rows = TRUE)

# check if microbiome::transform gives us the same
#pseq.comp <- microbiome::transform(biom.16s, "clr")
#d.clr.melt <- melt(otu_table(pseq.comp), id.vars = colnames(otu_table(pseq.comp))) %>% dplyr::rename("Sample" = "Var2", "clr" = "value")

# visualize clr distribution
d.clr.melt <- melt(otu_table(biom.16s.clr), id.vars = colnames(otu_table(biom.16s.clr))) %>% dplyr::rename("Sample" = "Var2", "clr" = "value")
d.clr.melt <- inner_join(d.clr.melt, meta, by = "Sample")
pdf(paste0(wdir,"fig/clr-distr-16s.pdf"), width = 14, height = 10)
ggplot(d.clr.melt, aes(clr, Sample, fill=Inoculum)) + geom_boxplot(outliers = F) + coord_flip() + facet_wrap(~ Culture_setup, ncol = 2, scales = "free_x") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
dev.off()

Most abundant genera

Maximal relative abundances (>3%) for each Culture setup aggregated for each genus

  temp.16s <- merge_samples(biom.16s, "Culture_setup") %>% tax_glom(taxrank = "Genus") %>% transform_sample_counts(function(x) {x / sum(x)}) %>% psmelt() %>% filter(Abundance>0.03) %>% arrange(Genus) %>% dplyr::arrange(desc(Abundance))
  knitr::kable(as.matrix(temp.16s[,c(2,3,13:17)]))

Sample	Abundance	Phylum	Class	Order	Family	Genus
light+clean_gas	0.33453718	Bacteroidota	Bacteroidia	Chitinophagales	Saprospiraceae	Aureispira
dark+clean_gas	0.32764094	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Marivita
ocean	0.27227980	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Lentibacter
light+used_gas	0.26110341	Bacteroidota	Bacteroidia	Chitinophagales	Saprospiraceae	Aureispira
light+used_gas	0.21513669	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Marivita
light+clean_gas	0.14536840	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Marivita
light+used_gas	0.13484422	Cyanobacteria	Cyanobacteriia	Cyanobacteriales	Geitlerinemaceae	Geitlerinema_PCC-7105
dark+clean_gas	0.12307436	Proteobacteria	Gammaproteobacteria	Nitrosococcales	Methylophagaceae	Methylophaga
light+clean_gas	0.10873896	Cyanobacteria	Cyanobacteriia	Cyanobacteriales	Geitlerinemaceae	Geitlerinema_PCC-7105
ocean	0.09912371	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	NS3a_marine_group
dark+clean_gas	0.08184165	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylomicrobium
light+used_gas	0.07883621	Bacteroidota	Bacteroidia	Chitinophagales	Saprospiraceae	Saprospiraceae_uncultured
ocean	0.05321352	Bacteroidota	Bacteroidia	Flavobacteriales	Cryomorphaceae	Cryomorphaceae_uncultured
ocean	0.05241535	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Planktomarina
light+clean_gas	0.04592369	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Flavobacterium
light+clean_gas	0.04436378	Bacteroidota	Bacteroidia	Chitinophagales	Saprospiraceae	Saprospiraceae_uncultured
light+used_gas	0.04364331	Bacteroidota	Bacteroidia	Chitinophagales	Saprospiraceae	Lewinella
light+clean_gas	0.03774251	Bacteroidota	Bacteroidia	Chitinophagales	Saprospiraceae	Lewinella
light+used_gas	0.03367726	Proteobacteria	Alphaproteobacteria	Rhodospirillales	Thalassospiraceae	Thalassospira
light+clean_gas	0.03128879	Proteobacteria	Alphaproteobacteria	Rhodospirillales	Thalassospiraceae	Thalassospira

Community typing with Dirichlet Multinomial Mixtures

set.seed(12345)

#https://microbiome.github.io/tutorials/DMM.html
library(DirichletMultinomial)

# prefilter relative abundances for prevalent taxa
biom.16s.rel.genus <- biom.16s.rel %>% tax_glom(taxrank = "Genus")
biom.16s.genus <- biom.16s %>% tax_glom(taxrank = "Genus")
taxa <- core_members(biom.16s.rel.genus, detection = 1/200, prevalence = 6/30) #table(meta$Culture_setup) -> minimum no. samples of dark-light-ocean: 6
pseq <- prune_taxa(taxa, biom.16s.genus)
dat <- abundances(pseq)
#rownames(dat) <- paste(tax_table(pseq)[,"Genus"], tax_table(pseq)[,"Species"], sep = "::")
rownames(dat) <- tax_table(pseq)[,"Genus"]
count <- as.matrix(t(dat))

# fit the DMM model
fit <- lapply(1:7, dmn, count = count, verbose=TRUE, seed = 1234567)

## dmn, k=1

##   Soft kmeans
##   Expectation Maximization setup
##   Expectation Maximization
##   Hessian

## dmn, k=2

##   Soft kmeans
##     iteration 10 change 0.000000
##   Expectation Maximization setup
##   Expectation Maximization
##   Hessian

## dmn, k=3

##   Soft kmeans
##     iteration 10 change 0.000007
##   Expectation Maximization setup
##   Expectation Maximization
##   Hessian

## dmn, k=4

##   Soft kmeans
##     iteration 10 change 0.012249
##   Expectation Maximization setup
##   Expectation Maximization
##   Hessian

## dmn, k=5

##   Soft kmeans
##   Expectation Maximization setup
##   Expectation Maximization
##   Hessian

## dmn, k=6

##   Soft kmeans
##     iteration 10 change 0.000001
##   Expectation Maximization setup
##   Expectation Maximization
##   Hessian

## dmn, k=7

##   Soft kmeans
##   Expectation Maximization setup
##   Expectation Maximization
##   Hessian

# Check model fit with different number of mixture components using standard information criteria
lplc <- base::sapply(fit, DirichletMultinomial::laplace) # AIC / BIC / Laplace
aic  <- base::sapply(fit, DirichletMultinomial::AIC) # AIC / BIC / Laplace
bic  <- base::sapply(fit, DirichletMultinomial::BIC) # AIC / BIC / Laplace
plot(lplc, type="b", xlab="Number of Dirichlet Components", ylab="Model Fit", ylim=c(min(c(lplc,aic,bic), na.rm = T), max(c(lplc,aic,bic), na.rm = T))); lines(aic, type="b", lty = 2); lines(bic, type="b", lty = 3)

# pick the optimal model
best <- fit[[which.min(unlist(lplc))]]

# mixture parameters pi and theta
mixturewt(best)

##          pi    theta
## 1 0.3750003 19.58217
## 2 0.3749997 18.62693
## 3 0.2500000 65.94614

# sample-component assignments
ass <- apply(mixture(best), 1, which.max)
temp.meta <- meta
temp.meta$cluster <- ass
knitr::kable(as.matrix(temp.meta %>% dplyr::select(Culture_setup, Inoculum, Time_point, cluster)))

Culture_setup	Inoculum	Time_point	cluster
ocean	Helsingor	0	3
light+clean_gas	Helsingor	7	1
light+clean_gas	Helsingor	14	2
ocean	Helsingor	0	3
ocean	Hornbaek-deep	0	3
ocean	Hornbaek-deep	0	3
light+clean_gas	Hornbaek-deep	5	1
light+clean_gas	Hornbaek-deep	7	1
light+clean_gas	Hornbaek-deep	11	1
light+clean_gas	Hornbaek-deep	11	1
light+used_gas	Hornbaek-deep	5	1
light+used_gas	Hornbaek-deep	7	1
light+used_gas	Hornbaek-deep	11	1
light+used_gas	Hornbaek-deep	11	1
dark+clean_gas	Hornbaek-deep	5	2
dark+clean_gas	Hornbaek-deep	7	2
dark+clean_gas	Hornbaek-deep	11	2
dark+clean_gas	Hornbaek-deep	11	2
ocean	Hornbaek-shallow	0	3
ocean	Hornbaek-shallow	0	3
dark+clean_gas	Hornbaek-shallow	5	2
dark+clean_gas	Hornbaek-shallow	7	2
dark+clean_gas	Hornbaek-shallow	11	2
dark+clean_gas	Hornbaek-shallow	11	2

# Contribution of each taxonomic group to each component
for (k in seq(ncol(fitted(best)))) {
  d <- melt(fitted(best))
  colnames(d) <- c("OTU", "cluster", "value")
  d <- subset(d, cluster == k) %>%
     # Arrange OTUs by assignment strength
     arrange(value) %>%
     mutate(OTU = factor(OTU, levels = unique(OTU))) %>%
     # Only show the most important drivers
     filter(abs(value) > quantile(abs(value), 0.8))     

  p <- ggplot(d, aes(x = OTU, y = value)) +
       geom_bar(stat = "identity") +
       coord_flip() +
       labs(title = paste("Top drivers: community type", k))
  print(p)
}

Cyanobacteria and Methanotrophs

Subset cyanobacteria and methanotrophs

#get cyanobacteria
biom.cyano <- subset_taxa(biom.16s, Phylum=="Cyanobacteria")
biom.cyano.clr <- subset_taxa(biom.16s.clr, Phylum=="Cyanobacteria")
biom.cyano.rel <- subset_taxa(biom.16s.rel, Phylum=="Cyanobacteria")
biom.cyano.rel

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 53 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 53 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 53 tips and 52 internal nodes ]

tax_stats(biom.cyano.rel)

## [1] "biom.cyano.rel: No. reads 1.7280778219388, No. species 53, No. genus 19, No. family 12"

#get methanotrophs
biom.mob <- subset_taxa(biom.16s, Phylum=="Methanotroph")
biom.mob.clr <- subset_taxa(biom.16s.clr, Phylum=="Methanotroph")
biom.mob.rel <- subset_taxa(biom.16s.rel, Phylum=="Methanotroph")
biom.mob.rel

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 22 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 22 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 22 tips and 21 internal nodes ]

tax_stats(biom.mob.rel)

## [1] "biom.mob.rel: No. reads 1.34902775524237, No. species 22, No. genus 8, No. family 1"

Cyanobacteria

cyanobac <- tax_glom(biom.cyano.rel, taxrank = "Phylum")
p1 <- plot_bar(cyanobac) + geom_bar(stat = "identity", position = "stack", fill = "dark green") + ggtitle("Cyanobacteria relative abundance") + theme(axis.text.x =element_text(angle=90,hjust=1)) + scale_x_discrete(limits = sampleid.16s)
cyanobac <- tax_glom(biom.cyano.clr, taxrank = "Phylum")
p2 <- plot_bar(cyanobac) + geom_bar(stat = "identity", position = "stack", fill = "dark green") + ggtitle("Cyanobacteria centered log-ratio") + theme(axis.text.x =element_text(angle=90,hjust=1)) + scale_x_discrete(limits = sampleid.16s)
plot(ggarrange(p1, p2, ncol = 2, nrow = 1))

#relative abundance of all bacteria
taxa.cyano <- biom.cyano.rel %>% tax_glom(taxrank = "Family") %>%
                                psmelt() %>%
                                filter(Abundance>0.001) %>%
                                arrange(Family)
no.taxa <- length(levels(as.factor(taxa.cyano$Family)))
getPalette = colorRampPalette(brewer.pal(8, "Set1"))
my.palette = getPalette(no.taxa)
#taxa.cyano$Family <- fct_reorder2(taxa.cyano$Family, taxa.cyano$Abundance, taxa.cyano$Time_point)
taxa.cyano$Culture_setup <- factor(taxa.cyano$Culture_setup, levels = c("ocean", "light+clean_gas", "light+used_gas", "dark+clean_gas"))
taxa.cyano$Time_point <- factor(taxa.cyano$Time_point, levels = c("0","5","7","11","14"))
pdf(paste0(wdir, "fig/family-over-time-cyano-rel2bact.pdf"), width = 13, height = 13)
ggplot(taxa.cyano, aes(x = Time_point, y = Abundance, fill = Family)) +  geom_bar(stat = "identity") + facet_grid(vars(Inoculum),vars(Culture_setup), scales = "free_x", space = "free_x") + scale_fill_manual(values = my.palette) + ylab("Relative abundance")
dev.off()

#relative abundance of all cyanobacteria
#set counts to zero if relative abundance is less than 0.001
filterfun1 = function(x){
  x[(x / sum(x)) < (0.0001)] <- 0
  return(x)
}
biom.16s.sub <- transform_sample_counts(biom.16s, fun = filterfun1)
biom.cyano.sub <- subset_taxa(biom.16s.sub, Phylum=="Cyanobacteria")
#samples with the same label should be merged for plotting of relative abundances for Culture_setup-Inoculum-Time_point categories
biom.cyano.nr <- biom.cyano.sub
sample_data(biom.cyano.nr) <- sample_data(biom.cyano.nr)[,c("Label", "Culture_setup", "Inoculum", "Time_point")]
biom.cyano.nr <- merge_samples(biom.cyano.nr, "Label") #merge replicates
temp.cyano.nr.sample_data <- sample_data(biom.cyano)[,c("Label", "Culture_setup", "Inoculum", "Time_point")]
#https://rdrr.io/github/kstagaman/phyloseqCompanion/src/R/sample_data_frame.R
sample.data.frame <- function(ps) {
  return(as(phyloseq::sample_data(ps), "data.frame"))
}
temp.cyano.nr.sample_data <- sample.data.frame(temp.cyano.nr.sample_data) %>% distinct()
rownames(temp.cyano.nr.sample_data) <- temp.cyano.nr.sample_data$Label
sample_data(biom.cyano.nr) <- temp.cyano.nr.sample_data
taxa.cyano <- biom.cyano.nr %>% tax_glom(taxrank = "Family") %>%
                                transform_sample_counts(function(x) {x / sum(x)}) %>%
                                psmelt() %>%
                                filter(Abundance>0.01) %>%
                                arrange(Family)
no.taxa <- length(levels(as.factor(taxa.cyano$Family)))
getPalette = colorRampPalette(brewer.pal(8, "Set1"))
my.palette = getPalette(no.taxa)
#taxa.cyano$Family <- fct_reorder2(taxa.cyano$Family, taxa.cyano$Abundance, taxa.cyano$Time_point)
taxa.cyano$Culture_setup <- factor(taxa.cyano$Culture_setup, levels = c("ocean", "light+clean_gas", "light+used_gas", "dark+clean_gas"))
taxa.cyano$Time_point <- factor(taxa.cyano$Time_point, levels = c("0","5","7","11","14"))
pdf(paste0(wdir, "fig/family-over-time-cyano-rel2cyano.pdf"), width = 13, height = 13)
ggplot(taxa.cyano, aes(x = Time_point, y = Abundance, fill = Family)) +  geom_bar(stat = "identity") + facet_grid(vars(Inoculum),vars(Culture_setup), scales = "free_x", space = "free_x") + scale_fill_manual(values = my.palette) + ylab("Relative abundance")
dev.off()

`

Methanotroph

mob <- tax_glom(biom.mob.rel, taxrank = "Phylum")
p1 <- plot_bar(mob) + geom_bar(stat = "identity", position = "stack", fill = "dark green") + ggtitle("Methanotroph relative abundance") + theme(axis.text.x =element_text(angle=90,hjust=1)) + scale_x_discrete(limits = sampleid.16s)
mob <- tax_glom(biom.mob.clr, taxrank = "Phylum")
p2 <- plot_bar(mob) + geom_bar(stat = "identity", position = "stack", fill = "dark green") + ggtitle("Methanotroph centered log-ratio") + theme(axis.text.x =element_text(angle=90,hjust=1)) + scale_x_discrete(limits = sampleid.16s)
plot(ggarrange(p1, p2, ncol = 2, nrow = 1))

#relative abundance of all bacteria
taxa.mob <- biom.mob.rel %>% tax_glom(taxrank = "Genus") %>%
                                psmelt() %>%
                                filter(Abundance>0.0001) %>%
                                arrange(Genus)
no.taxa <- length(levels(as.factor(taxa.mob$Genus)))
getPalette = colorRampPalette(brewer.pal(8, "Set2"))
my.palette = getPalette(no.taxa)
#taxa.mob$Family <- fct_reorder2(taxa.mob$Family, taxa.mob$Abundance, taxa.mob$Time_point)
taxa.mob$Culture_setup <- factor(taxa.mob$Culture_setup, levels = c("ocean", "light+clean_gas", "light+used_gas", "dark+clean_gas"))
taxa.mob$Time_point <- factor(taxa.mob$Time_point, levels = c("0","5","7","11","14"))
pdf(paste0(wdir, "fig/genus-over-time-mob-rel2mob.pdf"), width = 13, height = 13)
ggplot(taxa.mob, aes(x = Time_point, y = Abundance, fill = Genus)) +  geom_bar(stat = "identity") + facet_grid(vars(Inoculum),vars(Culture_setup), scales = "free_x", space = "free_x") + scale_fill_manual(values = my.palette) + ylab("Relative abundance")
dev.off()

## png 
##   2

#relative abundance of all methanotrophs
#set counts to zero if relative abundance is less than 0.001
filterfun1 = function(x){
  x[(x / sum(x)) < (0.0001)] <- 0
  return(x)
}
biom.16s.sub <- transform_sample_counts(biom.16s, fun = filterfun1)
biom.mob.sub <- subset_taxa(biom.16s.sub, Phylum=="Methanotroph")
#samples with the same label should be merged for plotting of relative abundances for Culture_setup-Inoculum-Time_point categories
biom.mob.nr <- biom.mob.sub
sample_data(biom.mob.nr) <- sample_data(biom.mob.nr)[,c("Label", "Culture_setup", "Inoculum", "Time_point")]
biom.mob.nr <- merge_samples(biom.mob.nr, "Label") #merge replicates
temp.mob.nr.sample_data <- sample_data(biom.mob)[,c("Label", "Culture_setup", "Inoculum", "Time_point")]
#https://rdrr.io/github/kstagaman/phyloseqCompanion/src/R/sample_data_frame.R
sample.data.frame <- function(ps) {
  return(as(phyloseq::sample_data(ps), "data.frame"))
}
temp.mob.nr.sample_data <- sample.data.frame(temp.mob.nr.sample_data) %>% distinct()
rownames(temp.mob.nr.sample_data) <- temp.mob.nr.sample_data$Label
sample_data(biom.mob.nr) <- temp.mob.nr.sample_data
taxa.mob <- biom.mob.nr %>% tax_glom(taxrank = "Genus") %>%
                                transform_sample_counts(function(x) {x / sum(x)}) %>%
                                psmelt() %>%
                                filter(Abundance>0.01) %>%
                                arrange(Genus)
no.taxa <- length(levels(as.factor(taxa.mob$Genus)))
getPalette = colorRampPalette(brewer.pal(8, "Set2"))
my.palette = getPalette(no.taxa)
#taxa.mob$Family <- fct_reorder2(taxa.mob$Family, taxa.mob$Abundance, taxa.mob$Time_point)
taxa.mob$Culture_setup <- factor(taxa.mob$Culture_setup, levels = c("ocean", "light+clean_gas", "light+used_gas", "dark+clean_gas"))
taxa.mob$Time_point <- factor(taxa.mob$Time_point, levels = c("0","5","7","11","14"))
pdf(paste0(wdir, "fig/genus-over-time-mob-rel2all.pdf"), width = 13, height = 13)
ggplot(taxa.mob, aes(x = Time_point, y = Abundance, fill = Genus)) +  geom_bar(stat = "identity") + facet_grid(vars(Inoculum),vars(Culture_setup), scales = "free_x", space = "free_x") + scale_fill_manual(values = my.palette) + ylab("Relative abundance")
dev.off()

## png 
##   2

mob.light <- taxa.mob %>% filter(Culture_setup=="light+clean_gas" & Inoculum=="Hornbaek-deep" & Time_point==11) %>% group_by(Genus) %>% summarise(sum_relative_abundance = sum(Abundance)) %>% arrange(desc(sum_relative_abundance))
knitr::kable(as.matrix(mob.light))

Genus	sum_relative_abundance
Methylomicrobium	0.6426891
Methylomonas	0.3573109

Preprocess minimap2 abundancies

Different statistical methods (alpha and beta diversity, DEA) require raw read counts instead of probabilities as calculated by Emu.

#get minimap2 filtered raw counts
sample.biom <- paste0(devdir, "data/16s.minimap2.all_json_md.biom")
sample.tree <- paste0(devdir, "data/16s.minimap2.tree.nwk")
biom.m <- import_biom(sample.biom, treefilename = sample.tree)
colnames(tax_table(biom.m)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
tax.m <- data.frame(tax_table(biom.m))
tax.clean.m <- data.frame(row.names = row.names(tax.m),
                        Kingdom = str_replace(tax.m[,1], "d__", ""),
                        Phylum = str_replace(tax.m[,2], "p__", ""),
                        Class = str_replace(tax.m[,3], "c__", ""),
                        Order = str_replace(tax.m[,4], "o__", ""),
                        Family = str_replace(tax.m[,5], "f__", ""),
                        Genus = str_replace(tax.m[,6], "g__", ""),
                        Species = str_replace(tax.m[,7], "s__", ""),
                        stringsAsFactors = FALSE)
tax_table(biom.m) <- as.matrix(tax.clean.m)
biom.m <- fix_duplicate_tax(biom.m)

#filter for bacteria
biom.filter.silva.m <- subset_taxa(biom.m, Order!="Chloroplast")
biom.filter.silva.m <- subset_taxa(biom.filter.silva.m, Family!="Mitochondria")
biom.16s.m <- subset_taxa(biom.filter.silva.m, Kingdom=="Bacteria")
biom.16s.m <- prune_samples(!(sample_names(biom.16s.m)) %in% c("flaregas-4","flaregas-5","flaregas-6","flaregas-7","flaregas-8","flaregas-9"), biom.16s.m)

#replace phylum to methanotrophs
tax <- data.frame(tax_table(biom.16s.m))
tax.clean <- tax %>% mutate(Phylum = if_else(Genus %in% methanotrophs.genus, "Methanotroph", Phylum))
tax.clean <- tax.clean %>% mutate(Phylum = if_else(Family %in% methanotrophs.family, "Methanotroph", Phylum))
tax_table(biom.16s.m) <- as.matrix(tax.clean)

#preprocess phyloseq object
biom.16s.m.unfiltered <- biom.16s.m
prev0.m = apply(X = otu_table(biom.16s.m),
              MARGIN = ifelse(taxa_are_rows(biom.16s.m), yes = 1, no = 2),
              FUN = function(x){sum(x > 0)})
prevdf.m = data.frame(Prevalence = prev0.m, TotalAbundance = taxa_sums(biom.16s.m), tax_table(biom.16s.m))
ggplot(prevdf.m, aes(Prevalence)) + geom_bar()

prevalenceThreshold.m = 0.05 * nsamples(biom.16s.m)
biom.16s.m = prune_taxa((prev0.m > prevalenceThreshold.m), biom.16s.m)
biom.16s.m

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 9641 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 9641 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 9641 tips and 9639 internal nodes ]

#tax_stats(biom.16s.m)

#centered log-ratio transformation
d.n0 <- cmultRepl(otu_table(biom.16s.m), method="CZM", label=0, z.warning=0.99)

## No. adjusted imputations:  182271

dim(d.n0)

## [1] 9641   24

d.clr <- apply(d.n0, 2, function(x) log(x) - mean(log(x)))
biom.16s.m.clr <- biom.16s.m
otu_table(biom.16s.m.clr) <- otu_table(d.clr, taxa_are_rows = TRUE)

#transform to relative abundance
biom.16s.m.rel <- biom.16s.m %>% transform_sample_counts(function(x) {x / sum(x)})

#genus
#get genus counts
biom.16s.m.genus <- biom.16s.m %>% tax_glom(taxrank = "Genus")

#filter genera that exist also running Emu on minimap2 counts
biom.16s.genus.minimap2emu <- tax_table(biom.16s.m.genus) %>% as.data.frame() %>% dplyr::select(Genus) %>% rownames_to_column(var = "SpeciesID") %>% inner_join(tax_table(biom.16s.genus) %>% as.data.frame() %>% dplyr::select(Genus), by = "Genus")

Alpha diversity measure

This should not be calculated with low abundancy removed counts (so here we should use minimap2 abundancies).

• Observed .. counts up the number of different taxa you observe in a sample at a given taxonomic level
• Chao0 .. species richness (similar to observed)
• Shannon .. estimator for both species richness and evenness; higher means higher uncertainty of predicting which species you would see next if you were to look at another read from this sample
• Simpson .. Gini-Simpson Index; similar idea to Shannon; based on the probability that two entities (microbes, or reads) taken from the sample at random are of different types

Sample richness in the different culture setups and different time points (emu)

#https://docs.onecodex.com/en/articles/4136553-alpha-diversity
alpha_meas = c("Observed", "Chao1", "Shannon", "Simpson")

p <- plot_richness(biom.16s.int, "Culture_setup", "Time_point", measures=alpha_meas)
pdf(paste0(wdir, "fig/alpha-diversity-cultures-16s-emu.pdf"), width = 20)
p + geom_boxplot(data=p$data, aes(x=Culture_setup, y=value, color=NULL), alpha=0.1)
dev.off()

## png 
##   2

#boxplot_alpha(biom.16s.int, index = "chao1", x_var = "Culture_setup") + theme_minimal() + theme(legend.position = "none")

#https://microbiome.github.io/tutorials/Alphadiversity.html
tab <-microbiome::alpha(biom.16s.int, index = "all")

## Observed richness

## Other forms of richness

## Diversity

## Evenness

## Dominance

## Rarity

knitr::kable(as.matrix(tab))

	observed	chao1	diversity_inverse_simpson	diversity_gini_simpson	diversity_shannon	diversity_fisher	diversity_coverage	evenness_camargo	evenness_pielou	evenness_simpson	evenness_evar	evenness_bulla	dominance_dbp	dominance_dmn	dominance_absolute	dominance_relative	dominance_simpson	dominance_core_abundance	dominance_gini	rarity_log_modulo_skewness	rarity_low_abundance	rarity_rare_abundance
flaregas-1	364	364	23.253857	0.9569964	4.446192	54.875603	12	0.4496586	0.7539556	0.0638842	0.4972124	0.4198321	0.1641734	0.2377957	6837	0.1641734	0.0430036	0.2410614	0.9343243	2.061228	0.2384680	0.3259455
flaregas-2	205	205	18.365758	0.9455508	4.043574	31.805434	8	0.9390088	0.7596405	0.0895891	0.5112685	0.4449620	0.1812000	0.2720500	3624	0.1812000	0.0544492	0.1379000	0.9596532	2.061093	0.1390500	0.5427000
flaregas-3	71	71	25.615507	0.9609611	3.716464	11.734826	10	0.9940452	0.8718610	0.3607818	0.5741982	0.5767494	0.1161667	0.1860278	577	0.1161667	0.0390389	0.2158244	0.9792600	2.061423	0.0000000	0.6746527
flaregas-10	144	144	18.012854	0.9444841	3.950388	22.743772	10	0.8833613	0.7948765	0.1250893	0.5435716	0.4969869	0.1970528	0.2659508	2514	0.1970528	0.0555159	0.2497257	0.9669065	2.061357	0.0723468	0.2570936
flaregas-11	206	206	46.140687	0.9783272	4.620210	33.537379	22	0.9002330	0.8671767	0.2239839	0.5779056	0.5624584	0.0993126	0.1437657	1546	0.0993126	0.0216728	0.1593756	0.9436860	2.061391	0.1217319	0.3237618
flaregas-12	584	584	56.467394	0.9822907	5.240850	92.057685	32	0.6257597	0.8227522	0.0966907	0.5335078	0.4814331	0.0929103	0.1368791	4858	0.0929103	0.0177093	0.1365158	0.8728897	2.061315	0.3495324	0.4301069
flaregas-13	133	133	5.216084	0.8082853	2.844222	18.169120	3	0.9030886	0.5815989	0.0392187	0.3736430	0.3075903	0.4162777	0.4942386	11416	0.4162777	0.1917147	0.2092328	0.9850484	2.060952	0.0921091	0.1503792
flaregas-14	235	235	4.585795	0.7819353	2.753450	30.123017	2	0.7927532	0.5043331	0.0195140	0.3547464	0.2308487	0.4415997	0.5661598	32496	0.4415997	0.2180647	0.6596138	0.9820239	2.060846	0.1150883	0.1154280
flaregas-15	57	57	4.161387	0.7596955	2.295168	8.531925	2	1.0000000	0.5676821	0.0730068	0.4680488	0.3220544	0.4242380	0.6595494	2881	0.4242380	0.2403045	0.7754381	0.9919723	2.060550	0.0038286	0.0942424
flaregas-16	165	165	4.544816	0.7799691	2.788147	23.651645	2	0.7988789	0.5460589	0.0275443	0.4723405	0.3081889	0.4356846	0.5965224	11025	0.4356846	0.2200309	0.6885201	0.9805619	2.061274	0.1189093	0.1890535
flaregas-17	166	166	6.736979	0.8515655	2.973364	22.442065	3	0.8852945	0.5816453	0.0405842	0.3694369	0.2757042	0.3312828	0.4893099	12117	0.3312828	0.1484345	0.1825514	0.9831067	2.061229	0.1046861	0.1589567
flaregas-18	225	225	4.870536	0.7946838	3.019310	31.981029	2	0.8099476	0.5574694	0.0216468	0.4621243	0.3121526	0.4322940	0.5455873	15694	0.4322940	0.2053162	0.6230443	0.9738634	2.061262	0.1442541	0.1886018
flaregas-19	130	130	4.618498	0.7834794	2.620662	18.450632	2	0.6299565	0.5383962	0.0355269	0.4856388	0.3100489	0.3501701	0.6525232	7411	0.3501701	0.2165206	0.7419675	0.9840085	2.061261	0.0921376	0.1169911
flaregas-20	119	119	5.459180	0.8168223	2.861293	18.036618	2	0.6244773	0.5987067	0.0458755	0.5088171	0.3595004	0.3576836	0.5847843	4725	0.3576836	0.1831777	0.6996215	0.9824050	2.061229	0.0790310	0.1140045
flaregas-21	239	239	28.818345	0.9652999	4.357086	36.843558	14	0.9215666	0.7956022	0.1205789	0.4847310	0.4665360	0.1080745	0.2086128	2610	0.1080745	0.0347001	0.0801242	0.9508478	2.061303	0.1517184	0.2709731
flaregas-22	88	88	8.108700	0.8766757	3.195895	13.604011	4	0.9961611	0.7137938	0.0921443	0.4947832	0.4369549	0.3206944	0.4221106	2808	0.3206944	0.1233243	0.1146642	0.9829127	2.061244	0.0220420	0.1657149
flaregas-23	246	246	28.973900	0.9654862	4.460298	38.573434	16	0.8807671	0.8101780	0.1177801	0.5087523	0.4917949	0.1491681	0.1946247	3380	0.1491681	0.0345138	0.1389293	0.9456512	2.061378	0.1621872	0.3472792
flaregas-24	401	401	37.387162	0.9732529	4.798996	61.241357	20	0.7496266	0.8006385	0.0932348	0.4638321	0.4535202	0.1302554	0.1642840	5558	0.1302554	0.0267471	0.1326459	0.9216715	2.061103	0.2276307	0.4183033
flaregas-25	415	415	59.882167	0.9833005	5.184550	68.941248	36	0.6181893	0.8600383	0.1442944	0.6065978	0.5517742	0.0958120	0.1316840	2711	0.0958120	0.0166995	0.1515462	0.8895443	2.061407	0.3208694	0.4553101
flaregas-26	414	414	66.395387	0.9849387	5.217768	68.670111	38	0.6336633	0.8658951	0.1603753	0.6058342	0.5542743	0.0852522	0.1238530	2425	0.0852522	0.0150613	0.1263139	0.8883371	2.061407	0.3275444	0.4586746
flaregas-27	223	223	12.262297	0.9184492	3.740517	32.020001	6	0.8737199	0.6917696	0.0549879	0.4347468	0.3750955	0.2327347	0.3663969	7879	0.2327347	0.0815508	0.0862527	0.9665130	2.060918	0.1472204	0.2111420
flaregas-28	89	89	20.487573	0.9511899	3.778135	15.054468	10	0.9967909	0.8417111	0.2301974	0.5860519	0.5689563	0.1543455	0.2641543	856	0.1543455	0.0488101	0.2197981	0.9754881	2.061423	0.0059502	0.2068157
flaregas-29	216	216	2.195784	0.5445817	2.023634	27.457879	1	0.2482530	0.3764706	0.0101657	0.4024311	0.2137360	0.6720996	0.7063309	48123	0.6720996	0.4554183	0.7373780	0.9853563	2.060792	0.0992025	0.1068142
flaregas-30	116	116	2.523463	0.6037191	2.127191	15.683955	1	0.3663683	0.4474914	0.0217540	0.4277261	0.2617240	0.6249462	0.6705555	15963	0.6249462	0.3962809	0.6934581	0.9887809	2.061219	0.0790432	0.1013193

#show the dominant taxonomic group for each sample
temp <- sample_data(biom.16s.int)
temp$dominant <- tax_table(biom.16s.int)[dominant(biom.16s.int),"Species"]
knitr::kable(as.matrix(temp[,c("Culture_setup","dominant")]))

	Culture_setup	dominant
flaregas-1	ocean	Lentibacter_uncultured_Rhodobacteraceae
flaregas-2	light+clean_gas	Flavobacterium_uncultured_bacterium
flaregas-3	light+clean_gas	Marivita_Marivita_sp.
flaregas-10	ocean	Lentibacter_uncultured_Rhodobacteraceae
flaregas-11	ocean	Lentibacter_uncultured_Rhodobacteraceae
flaregas-12	ocean	Lentibacter_uncultured_Rhodobacteraceae
flaregas-13	light+clean_gas	Phormidiaceae_cyanobacterium
flaregas-14	light+clean_gas	Aureispira_marina
flaregas-15	light+clean_gas	Aureispira_marina
flaregas-16	light+clean_gas	Aureispira_marina
flaregas-17	light+used_gas	Phormidiaceae_cyanobacterium
flaregas-18	light+used_gas	Aureispira_marina
flaregas-19	light+used_gas	Marivita_Marivita_sp.
flaregas-20	light+used_gas	Marivita_Marivita_sp.
flaregas-21	dark+clean_gas	Methylophaga_thalassica
flaregas-22	dark+clean_gas	Methylomicrobium_uncultured_gamma
flaregas-23	dark+clean_gas	Methylophaga_thalassica
flaregas-24	dark+clean_gas	Methylophaga_thalassica
flaregas-25	ocean	Lentibacter_uncultured_Rhodobacteraceae
flaregas-26	ocean	Lentibacter_uncultured_Rhodobacteraceae
flaregas-27	dark+clean_gas	Methylomicrobium_uncultured_gamma
flaregas-28	dark+clean_gas	Saprospira_grandis
flaregas-29	dark+clean_gas	Marivita_Marivita_sp.
flaregas-30	dark+clean_gas	Marivita_Marivita_sp.

p <- plot_richness(biom.16s.int, "Time_point", "Culture_setup", measures=alpha_meas)
p$data$Time_point <- factor(p$data$Time_point, levels = c("0","5","7","11","14"))
pdf(paste0(wdir, "fig/alpha-diversity-times-16s-emu.pdf"), width = 20)
p + geom_boxplot(data=p$data, aes(x=Time_point, y=value, color=NULL), alpha=0.1)
dev.off()

## png 
##   2

#boxplot_alpha(biom.16s.int, index = "chao1", x_var = "Time_point") + theme_minimal() + theme(legend.position = "none")

Sample richness in the different culture setups and different time points (minimap2)

alpha_meas = c("Observed", "Chao1", "Shannon", "Simpson")

p <- plot_richness(biom.16s.m.unfiltered, "Culture_setup", "Time_point", measures=alpha_meas)
pdf(paste0(wdir, "fig/alpha-diversity-cultures-16s-minimap2.pdf"), width = 20)
p + geom_boxplot(data=p$data, aes(x=Culture_setup, y=value, color=NULL), alpha=0.1)
dev.off()

## png 
##   2

#boxplot_alpha(biom.16s.m.unfiltered, index = "chao1", x_var = "Culture_setup") + theme_minimal() + theme(legend.position = "none")

#https://microbiome.github.io/tutorials/Alphadiversity.html
tab <-microbiome::alpha(biom.16s.m.unfiltered, index = "all")

## Observed richness

## Other forms of richness

## Diversity

## Evenness

## Dominance

## Rarity

knitr::kable(as.matrix(tab))

	observed	chao1	diversity_inverse_simpson	diversity_gini_simpson	diversity_shannon	diversity_fisher	diversity_coverage	evenness_camargo	evenness_pielou	evenness_simpson	evenness_evar	evenness_bulla	dominance_dbp	dominance_dmn	dominance_absolute	dominance_relative	dominance_simpson	dominance_core_abundance	dominance_gini	rarity_log_modulo_skewness	rarity_low_abundance	rarity_rare_abundance
flaregas-1	4315	8445.548	15.286749	0.9345839	5.281947	1207.0809	20	0.8620797	0.6310681	0.0035427	0.4566371	0.2752696	0.2456530	0.2742429	10285	0.2456530	0.0654161	0.2467517	0.9617968	2.061421	0.3704261	0.4843795
flaregas-2	1987	4116.146	22.767845	0.9560784	4.924868	541.5041	15	0.9276232	0.6484883	0.0114584	0.4167156	0.2742059	0.1694604	0.2537076	3508	0.1694604	0.0439216	0.0453601	0.9827436	2.061421	0.3003720	0.6296314
flaregas-3	931	2259.231	31.718452	0.9684726	4.992567	314.0958	19	0.8662062	0.7303069	0.0340692	0.4691053	0.3637613	0.1391199	0.2011435	803	0.1391199	0.0315274	0.1476091	0.9875918	2.061421	0.3210326	0.7619543
flaregas-10	2298	5702.450	16.161425	0.9381243	5.209466	781.9143	21	0.9028564	0.6730754	0.0070328	0.5294223	0.3304788	0.2381905	0.2671336	3333	0.2381905	0.0618757	0.2400486	0.9718780	2.061423	0.3801901	0.4884585
flaregas-11	3563	8438.965	63.305943	0.9842037	6.300107	1332.8309	72	0.9212405	0.7703388	0.0177676	0.5320364	0.3890814	0.1070872	0.1415777	1925	0.1070872	0.0157963	0.1231086	0.9461639	2.061423	0.5018914	0.5826658
flaregas-12	6587	12834.554	55.341901	0.9819305	6.349247	1948.3011	68	0.8285402	0.7220918	0.0084017	0.4451133	0.3134241	0.1155695	0.1570329	6394	0.1155695	0.0180695	0.1242635	0.9302238	2.061421	0.4982829	0.5931750
flaregas-13	1349	3023.204	15.658169	0.9361356	3.979428	296.1536	6	0.9714376	0.5521524	0.0116072	0.3485270	0.1775566	0.1786022	0.3019877	4978	0.1786022	0.0638644	0.1423292	0.9935643	2.061400	0.1590126	0.2386625
flaregas-14	2105	4383.004	7.581352	0.8680974	3.590428	403.2499	4	0.8841795	0.4692099	0.0036016	0.3184866	0.1423329	0.3376338	0.4283211	25045	0.3376338	0.1319026	0.1830327	0.9933766	2.061407	0.1515409	0.1941546
flaregas-15	774	1989.101	8.703364	0.8851019	3.715408	218.5239	4	0.9642672	0.5585760	0.0112447	0.4723425	0.2480729	0.3041758	0.3906932	2229	0.3041758	0.1148981	0.2463155	0.9936617	2.061412	0.2083788	0.2194323
flaregas-16	1524	3329.817	7.788971	0.8716133	3.788956	353.0381	4	0.9236048	0.5169748	0.0051109	0.3860155	0.2056010	0.3329758	0.4226070	8693	0.3329758	0.1283867	0.1938944	0.9915259	2.061411	0.2021680	0.2710384
flaregas-17	1745	4175.615	15.314810	0.9347037	3.990306	379.7383	6	0.9748314	0.5345704	0.0087764	0.3631715	0.1668146	0.1670339	0.3014777	6217	0.1670339	0.0652963	0.1302526	0.9923364	2.061407	0.1630306	0.2489790
flaregas-18	1756	3646.043	7.730096	0.8706355	3.864414	384.5041	4	0.8923293	0.5172696	0.0044021	0.3496265	0.1959522	0.3388035	0.4228217	12408	0.3388035	0.1293645	0.1602818	0.9912086	2.061416	0.1974988	0.2674276
flaregas-19	1417	3100.050	11.228064	0.9109375	3.857661	337.9353	4	0.9173522	0.5316295	0.0079238	0.4020619	0.2038845	0.2263304	0.3497255	4989	0.2263304	0.0890625	0.3596607	0.9918641	2.061417	0.2070045	0.2215216
flaregas-20	1362	2793.256	16.655080	0.9399583	4.344194	370.0748	6	0.9047069	0.6019633	0.0122284	0.4394832	0.2576154	0.1592927	0.2830782	2279	0.1592927	0.0600417	0.3619207	0.9889622	2.061418	0.2472216	0.2664430
flaregas-21	2055	3721.157	50.490371	0.9801942	5.217657	529.3173	21	0.8893018	0.6840110	0.0245695	0.3697369	0.2648962	0.0588563	0.1155665	1481	0.0588563	0.0198058	0.0583794	0.9831396	2.061410	0.3097008	0.4206176
flaregas-22	1066	2646.504	15.908887	0.9371421	4.407192	308.5543	7	0.9403829	0.6321574	0.0149239	0.4401330	0.2874822	0.1723409	0.3289279	1630	0.1723409	0.0628579	0.0704166	0.9898860	2.061421	0.2544936	0.3418270
flaregas-23	2050	4228.650	53.969719	0.9814711	5.391625	535.9633	26	0.9035793	0.7070432	0.0263267	0.3618201	0.2839977	0.0826605	0.1485890	1986	0.0826605	0.0185289	0.0701740	0.9815394	2.061410	0.3318072	0.5316324
flaregas-24	2739	4807.189	59.213644	0.9831120	5.514769	649.5860	31	0.8502754	0.6967185	0.0216187	0.3258944	0.2588662	0.0825517	0.1426794	3582	0.0825517	0.0168880	0.0679634	0.9790365	2.061403	0.3179000	0.5762024
flaregas-25	5511	11670.811	72.584285	0.9862229	6.672463	1934.4862	114	0.8167415	0.7745617	0.0131708	0.5002829	0.3782221	0.1042931	0.1304458	3282	0.1042931	0.0137771	0.1259017	0.9215018	2.061423	0.5673202	0.6504496
flaregas-26	5460	10824.639	89.489842	0.9888255	6.720246	1912.8569	117	0.8065763	0.7809513	0.0163901	0.4993797	0.3821001	0.0886553	0.1205073	2775	0.0886553	0.0111745	0.1097409	0.9214100	2.061423	0.5704610	0.6566883
flaregas-27	2216	4300.333	21.730904	0.9539826	4.530733	528.6524	8	0.7515981	0.5881427	0.0098064	0.3909616	0.2174082	0.1296861	0.2391033	4466	0.1296861	0.0460174	0.0645527	0.9863596	2.061412	0.2435752	0.3862996
flaregas-28	1019	2374.803	4.567576	0.7810655	2.909577	236.4565	2	0.9764584	0.4200599	0.0044824	0.4329404	0.1798135	0.3684029	0.6516478	6394	0.3684029	0.2189345	0.0467274	0.9949559	2.061419	0.1496889	0.1669740
flaregas-29	1929	4190.534	8.079105	0.8762239	3.281015	364.7623	3	0.9760829	0.4337238	0.0041882	0.3402545	0.1270341	0.2398770	0.4579773	17239	0.2398770	0.1237761	0.6882531	0.9946194	2.061414	0.1338463	0.1716250
flaregas-30	1140	2177.682	9.508020	0.8948256	3.495112	243.9462	4	0.9840635	0.4965506	0.0083404	0.3623653	0.1691140	0.2259636	0.4202265	5845	0.2259636	0.1051744	0.6286388	0.9949990	2.061408	0.1481811	0.1886574

#show the dominant taxonomic group for each sample
temp <- sample_data(biom.16s.m.unfiltered)
temp$dominant <- tax_table(biom.16s.m.unfiltered)[dominant(biom.16s.m.unfiltered),"Species"]
knitr::kable(as.matrix(temp[,c("Culture_setup","dominant")]))

	Culture_setup	dominant
flaregas-1	ocean	Lentibacter_algarum
flaregas-2	light+clean_gas	Flavobacterium_uncultured_bacterium
flaregas-3	light+clean_gas	Psychroflexus_uncultured_organism
flaregas-10	ocean	Lentibacter_algarum
flaregas-11	ocean	Lentibacter_algarum
flaregas-12	ocean	Lentibacter_algarum
flaregas-13	light+clean_gas	Phormidiaceae_cyanobacterium
flaregas-14	light+clean_gas	Aureispira_marina
flaregas-15	light+clean_gas	Aureispira_marina
flaregas-16	light+clean_gas	Aureispira_marina
flaregas-17	light+used_gas	Saprospiraceae_uncultured_uncultured_marine
flaregas-18	light+used_gas	Aureispira_marina
flaregas-19	light+used_gas	Aureispira_marina
flaregas-20	light+used_gas	Aureispira_marina
flaregas-21	dark+clean_gas	NS3a_marine_group_uncultured_bacterium
flaregas-22	dark+clean_gas	Methylomicrobium_uncultured_bacterium
flaregas-23	dark+clean_gas	Methylophaga_thalassica
flaregas-24	dark+clean_gas	Methylophaga_thalassica
flaregas-25	ocean	Lentibacter_algarum
flaregas-26	ocean	Lentibacter_algarum
flaregas-27	dark+clean_gas	Methylomicrobium_uncultured_bacterium
flaregas-28	dark+clean_gas	Methylomicrobium_uncultured_bacterium
flaregas-29	dark+clean_gas	Marivita_Marivita_sp.
flaregas-30	dark+clean_gas	Marivita_Marivita_sp.

p <- plot_richness(biom.16s.m.unfiltered, "Time_point", "Culture_setup", measures=alpha_meas)
p$data$Time_point <- factor(p$data$Time_point, levels = c("0","5","7","11","14"))
pdf(paste0(wdir, "fig/alpha-diversity-times-16s-minimap2.pdf"), width = 20)
p + geom_boxplot(data=p$data, aes(x=Time_point, y=value, color=NULL), alpha=0.1)
dev.off()

## png 
##   2

#boxplot_alpha(biom.16s.m.unfiltered, index = "chao1", x_var = "Time_point") + theme_minimal() + theme(legend.position = "none")

Beta diversity measure

This should not be calculated with low abundancy removed counts (so here we should use minimap2 abundancies).

PCA (emu)

# apply a singular value decomposition to the dataset
# do not use princomp function in R!!
pcx <- prcomp(t(otu_table(biom.16s.clr)), scale. = T)

# generate a scree plot
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/
fviz_eig(pcx, addlabels = TRUE)

# generate biplot (check https://colorbrewer2.org/ for selecting colors)
pdf(paste0(wdir,"fig/clr-pca-16s-emu.pdf"), width = 30, height = 8)
p1 <- ggbiplot(pcx, groups = meta$Culture_setup, var.axes = F, choices = c(1,3), obs.scale = 1, var.scale = 1, ellipse = T, circle = T, ellipse.fill = F, labels = meta$Label, labels.size = 5) + scale_color_brewer(type = "div", palette = 1) + labs(color = "Group")
p2 <- ggbiplot(pcx, groups = meta$Time_point, var.axes = F, choices = c(1,4), obs.scale = 1, var.scale = 1, ellipse = T, circle = T, ellipse.fill = F, labels = meta$Label, labels.size = 5) + scale_color_manual(values=c('#e41a1c','#377eb8','#4daf4a','#984ea3','#ff7f00')) + labs(color = "Group")
plot(ggarrange(p1, p2, ncol = 2, nrow = 1))

## Too few points to calculate an ellipse

dev.off()

PCA (minimap2)

#biplot
pcx <- prcomp(t(otu_table(biom.16s.m.clr)), scale. = T)
fviz_eig(pcx, addlabels = TRUE)

# generate biplot (check https://colorbrewer2.org/ for selecting colors)
p1 <- ggbiplot(pcx, groups = meta$Culture_setup, var.axes = F, choices = c(1,2), obs.scale = 1, var.scale = 1, ellipse = T, circle = T, ellipse.fill = F, labels = meta$Label, labels.size = 5) + scale_color_brewer(type = "div", palette = 1) + labs(color = "Group")
p2 <- ggbiplot(pcx, groups = meta$Time_point, var.axes = F, choices = c(1,2), obs.scale = 1, var.scale = 1, ellipse = T, circle = T, ellipse.fill = F, labels = meta$Label, labels.size = 5) + scale_color_manual(values=c('#e41a1c','#377eb8','#4daf4a','#984ea3','#ff7f00')) + labs(color = "Group")
pdf(paste0(wdir,"fig/clr-pca-16s-minimap2.pdf"), width = 30, height = 8)
plot(ggarrange(p1, p2, ncol = 2, nrow = 1))

## Too few points to calculate an ellipse

dev.off()

Hierarchical Clustering

# average-link clustering
# Use unifrac on original compositional data: If your data is truly compositional (relative abundances sum to 1), consider using unifrac directly on the untransformed data
biom.hclust <- hclust(phyloseq::distance(biom.16s.m, method = "unifrac"), method="average")
biom.hclust$labels <- sample_data(biom.16s)$Label
plot(biom.hclust)

#biom.hclust <- hclust(phyloseq::distance(biom.16s, "wunifrac"), method="average")
#biom.hclust$labels <- sample_data(biom.16s)$Label
#plot(biom.hclust)
#biom.hclust <- hclust(phyloseq::distance(biom.16s, "bray"), method="average")
#biom.hclust$labels <- sample_data(biom.16s)$Label
#plot(biom.hclust)

Multivariate community comparison

PERMANOVA

This should not be calculated with low abundancy removed counts (so here we should use minimap2 abundancies).

#https://microbiome.github.io/tutorials/PERMANOVA.html
plot_landscape(biom.16s.m.rel, method = "NMDS", distance = "bray", col = "Culture_setup", size = 3)

plot_landscape(biom.16s.m.rel, method = "NMDS", distance = "bray", col = "Time_point", size = 3)

#https://deneflab.github.io/MicrobeMiseq/demos/mothur_2_phyloseq.html
set.seed(1)

# Calculate bray curtis distance matrix
#biom.16s.dist <- phyloseq::distance(biom.16s.m, method = "bray")
#biom.16s.dist <- phyloseq::distance(biom.16s.m, method = "unifrac")
biom.16s.dist <- phyloseq::distance(biom.16s.m, method = "wunifrac")

# make a data frame from the sample_data
#sampledf <- data.frame(sample_data(biom.16s.m))
sampledf <- meta(biom.16s.m)

#test the hypothesis that the different culture setups have different centroids
# Adonis test (H0 .. different cultures have the same centroid)
perm <- adonis2(biom.16s.dist ~ Culture_setup, data = sampledf)
perm

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = biom.16s.dist ~ Culture_setup, data = sampledf)
##               Df SumOfSqs      R2      F Pr(>F)    
## Culture_setup  3  0.43518 0.44116 5.2629  0.001 ***
## Residual      20  0.55125 0.55884                  
## Total         23  0.98643 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Homogeneity of dispersion test (H0 .. different cultures have the same dispersion)
beta <- betadisper(biom.16s.dist, sampledf$Culture_setup)
permutest(beta) #alternatively use anova(beta)

## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     3 0.021279 0.0070930 0.8965    999  0.441
## Residuals 20 0.158234 0.0079117

#test the hypothesis that the different time points have different centroids
# Adonis test (H0 .. different different time points have the same centroid)
perm <- adonis2(biom.16s.dist ~ Time_point, data = sampledf)
perm

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = biom.16s.dist ~ Time_point, data = sampledf)
##            Df SumOfSqs      R2     F Pr(>F)    
## Time_point  4  0.39174 0.39713 3.129  0.001 ***
## Residual   19  0.59468 0.60287                 
## Total      23  0.98643 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Homogeneity of dispersion test (H0 .. different time points have the same dispersion)
beta <- betadisper(biom.16s.dist, sampledf$Time_point)
permutest(beta)

## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
## Groups     4 0.054529 0.0136322 2.8825    999  0.042 *
## Residuals 19 0.089858 0.0047294                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Univariate community comparisons

ANOVA-like differential expression of genera (minimap2)

set.seed(12345)
# test differential expression: (light.clean_gas+light.used_gas) - dark.clean_gas
meta.sub.1 <- meta %>% dplyr::filter(Culture_setup == "light+clean_gas" | Culture_setup == "light+used_gas" | Culture_setup == "dark+clean_gas") %>% mutate(Culture_setup = replace(Culture_setup, Culture_setup=="light+used_gas", "light+clean_gas"))
mm <- model.matrix(~ Inoculum + Culture_setup, meta.sub.1)
counts.sub <- otu_table(biom.16s.m.genus)[,meta.sub.1$Sample]
x <- aldex.clr(counts.sub, mm, mc.samples=100, verbose=T, gamma=1e-3)

## integer matrix provided

## using all features for denominator

## operating in serial mode

## removed rows with sums equal to zero

## data format is OK

## dirichlet samples complete

## aldex.scaleSim: adjusting samples to reflect scale uncertainty.

## sampling from the default scale model.

glm.test <- aldex.glm(x, mm, fdr.method='BH')

## running tests for each MC instance:

## |------------(25%)----------(50%)----------(75%)----------|

glm.eff<- aldex.glm.effect(x)

## operating in serial mode

## sanity check complete

## rab.all  complete

## rab.win  complete

## rab of samples complete

## within sample difference calculated

## between group difference calculated

## group summaries calculated

## unpaired effect size calculated

## summarizing output

## operating in serial mode

## sanity check complete

## rab.all  complete

## rab.win  complete

## rab of samples complete

## within sample difference calculated

## between group difference calculated

## group summaries calculated

## unpaired effect size calculated

## summarizing output

## operating in serial mode

## sanity check complete

## rab.all  complete

## rab.win  complete

## rab of samples complete

## within sample difference calculated

## between group difference calculated

## group summaries calculated

## unpaired effect size calculated

## summarizing output

aldex.glm.plot(glm.test, eff=glm.eff, contrast='Culture_setuplight+clean_gas', type='MW', test='effect', cutoff.effect = 1)

aldex.glm.plot(glm.test, eff=glm.eff, contrast='Culture_setuplight+clean_gas', type='MA', test='effect', cutoff.effect = 1)

#aldex.glm.plot(glm.test, eff=glm.eff, contrast='Culture_setuplight+clean_gas', type='volcano', test='effect', cutoff.effect = 1)

#filter genera that where detected with Emu (filter low abundant genera)
glm.eff.lightVSdark <- glm.eff$`Culture_setuplight+clean_gas`[biom.16s.genus.minimap2emu$SpeciesID,]

# higher relative abundance under light
glm.eff.lightVSdark[glm.eff.lightVSdark$effect >= 1.5,] %>% nrow()

## [1] 12

knitr::kable(as.matrix(tax_table(biom.16s.m.genus)[glm.eff.lightVSdark[glm.eff.lightVSdark$effect >= 1.5,] %>% row.names(),]))

	Kingdom	Phylum	Class	Order	Family	Genus	Species
KY770139.1.1394	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodospirillales	Thalassospiraceae	Thalassospira	NA
EU672857.1.1492	Bacteria	Proteobacteria	Gammaproteobacteria	Cellvibrionales	Halieaceae	Congregibacter	NA
JN428535.1.1439	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Flammeovirgaceae	Flexithrix	NA
KU578816.1.1368	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	NS2b_marine_group	NA
AB245935.1.1459	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Saprospiraceae	Aureispira	NA
HM057705.1.1470	Bacteria	Cyanobacteria	Cyanobacteriia	Synechococcales	Cyanobiaceae	Cyanobium_PCC-6307	NA
KU951800.1.1271	Bacteria	Cyanobacteria	Cyanobacteriia	Synechococcales	Synechococcales_Incertae_Sedis	Schizothrix_LEGE_07164	NA
HM217045.1.1446	Bacteria	Cyanobacteria	Cyanobacteriia	Oxyphotobacteria_Incertae_Sedis	Oxyphotobacteria_Incertae_Sedis_Unknown_Family	Leptolyngbya_LEGE-06070	NA
EF654088.1.1482	Bacteria	Cyanobacteria	Cyanobacteriia	Cyanobacteriales	Geitlerinemaceae	Geitlerinema_PCC-7105	NA
KY379855.1.2514	Bacteria	Cyanobacteria	Cyanobacteriia	Cyanobacteriales	Microcystaceae	Synechocystis_PCC-6803	NA
KU951847.1.1466	Bacteria	Cyanobacteria	Cyanobacteriia	Cyanobacteriales	Cyanobacteriaceae	Symphothece_PCC-7002	NA
AB062106.1.1479	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Porphyrobacter	NA

# higher relative abundance in the dark
glm.eff.lightVSdark[glm.eff.lightVSdark$effect <= -1.5,] %>% nrow()

## [1] 16

knitr::kable(as.matrix(tax_table(biom.16s.m.genus)[glm.eff.lightVSdark[glm.eff.lightVSdark$effect <= -1.5,] %>% row.names(),]))

	Kingdom	Phylum	Class	Order	Family	Genus	Species
KC247325.1.1400	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Tropicibacter	NA
KY437089.1.1426	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Youngimonas	NA
AY546508.1.1498	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylomonadaceae_uncultured	NA
AY007296.1.1499	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylosarcina	NA
AB930603.1.1457	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylomicrobium	NA
CP002738.392605.394125	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylomonas	NA
KC963966.1.1493	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylobacter	NA
ATTE01000001.3195237.3196755	Bacteria	Proteobacteria	Gammaproteobacteria	Thiotrichales	Thiotrichaceae	Leucothrix	NA
AB681780.1.1463	Bacteria	Proteobacteria	Gammaproteobacteria	Nitrosococcales	Methylophagaceae	Methylophaga	NA
KC534174.1.1363	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Aequorivita	NA
JQ197367.1.1350	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Aquibacter	NA
LN875339.1.1446	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Flavobacteriaceae	NA
HQ144198.1.1471	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Mariniflexile	NA
FJ360684.1.1470	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Meridianimaribacter	NA
KJ475509.1.1443	Bacteria	Myxococcota	Polyangia	Nannocystales	Nannocystaceae	Nannocystaceae_uncultured	NA
AB607876.1.1417	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Antarctobacter	NA

temp.1 <- glm.eff.lightVSdark[abs(glm.eff.lightVSdark$effect) >= 1.5,] %>% rownames_to_column(var = "SpeciesID")
temp.2 <- tax_table(biom.16s.m.genus) %>% as.data.frame() %>% rownames_to_column(var = "SpeciesID")
res_sig.aldex <- inner_join(temp.1, temp.2, by = "SpeciesID")
pdf(paste0(wdir, "fig/aldex-light_clean_usedVSdark-scatter-genus-16s.pdf"), width = 35, height = 13)
ggplot(res_sig.aldex, aes(x=Genus, y=diff.btw, color=Phylum)) +
    geom_jitter(size=3, width = 0.2) +
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5)) + ylab("logFC")
dev.off()

## png 
##   2

for(i in glm.eff.lightVSdark[abs(glm.eff.lightVSdark$effect) >= 1,] %>% row.names()) {
  de.genus <- tax_table(biom.16s.m.genus)[i,"Genus"]
  p <- boxplot_abundance(biom.16s.m.genus, x = "Culture_setup", y = i) + scale_y_log10() + ggtitle(de.genus)
  print(p)
}

set.seed(12345)
# test differential expression: ((light.clean_gas+light.used_gas) - ocean) AND (dark.clean_gas - ocean)
meta.sub.1 <- meta %>% mutate(Culture_setup = replace(Culture_setup, Culture_setup=="light+used_gas", "light+clean_gas"))
meta.sub.1$Culture_setup <- factor(meta.sub.1$Culture_setup, levels = c("ocean", "dark+clean_gas", "light+clean_gas"))
mm <- model.matrix(~0 + Inoculum + Culture_setup, meta.sub.1)
counts.sub <- otu_table(biom.16s.m.genus)[,meta.sub.1$Sample]
x <- aldex.clr(counts.sub, mm, mc.samples=100, verbose=T, gamma=1e-3)

## integer matrix provided

## using all features for denominator

## operating in serial mode

## removed rows with sums equal to zero

## data format is OK

## dirichlet samples complete

## aldex.scaleSim: adjusting samples to reflect scale uncertainty.

## sampling from the default scale model.

glm.test.1 <- aldex.glm(x, mm, fdr.method='BH')

## running tests for each MC instance:

## |------------(25%)----------(50%)----------(75%)----------|

glm.eff.1<- aldex.glm.effect(x)

## operating in serial mode

## sanity check complete

## rab.all  complete

## rab.win  complete

## rab of samples complete

## within sample difference calculated

## between group difference calculated

## group summaries calculated

## unpaired effect size calculated

## summarizing output

## operating in serial mode

## sanity check complete

## rab.all  complete

## rab.win  complete

## rab of samples complete

## within sample difference calculated

## between group difference calculated

## group summaries calculated

## unpaired effect size calculated

## summarizing output

## operating in serial mode

## sanity check complete

## rab.all  complete

## rab.win  complete

## rab of samples complete

## within sample difference calculated

## between group difference calculated

## group summaries calculated

## unpaired effect size calculated

## summarizing output

## operating in serial mode

## sanity check complete

## rab.all  complete

## rab.win  complete

## rab of samples complete

## within sample difference calculated

## between group difference calculated

## group summaries calculated

## unpaired effect size calculated

## summarizing output

aldex.glm.plot(glm.test.1, eff=glm.eff.1, contrast='Culture_setuplight+clean_gas', type='MW', test='effect', cutoff.effect = 1.5)

aldex.glm.plot(glm.test.1, eff=glm.eff.1, contrast='Culture_setuplight+clean_gas', type='MA', test='effect', cutoff.effect = 1.5)

#aldex.glm.plot(glm.test.1, eff=glm.eff.1, contrast='Culture_setuplight+clean_gas', type='volcano', test='fdr', cutoff.pval = 0.1)
aldex.glm.plot(glm.test.1, eff=glm.eff.1, contrast='Culture_setupdark+clean_gas', type='MW', test='effect', cutoff.effect = 1.5)

aldex.glm.plot(glm.test.1, eff=glm.eff.1, contrast='Culture_setupdark+clean_gas', type='MA', test='effect', cutoff.effect = 1.5)

#aldex.glm.plot(glm.test.1, eff=glm.eff.1, contrast='Culture_setupdark+clean_gas', type='volcano', test='fdr', cutoff.pval = 0.1)

#filter genera that where detected with Emu (filter low abundant genera)
glm.eff.lightVSocean <- glm.eff.1$`Culture_setuplight+clean_gas`[biom.16s.genus.minimap2emu$SpeciesID,]
glm.eff.darkVSocean <- glm.eff.1$`Culture_setupdark+clean_gas`[biom.16s.genus.minimap2emu$SpeciesID,]

# higher relative abundance in light-incubated cultures than in the ocean
glm.eff.lightVSocean[glm.eff.lightVSocean$effect >= 1.5,] %>% nrow()

## [1] 9

knitr::kable(as.matrix(tax_table(biom.16s.m.genus)[glm.eff.lightVSocean[glm.eff.lightVSocean$effect >= 1.5,] %>% row.names(),]))

	Kingdom	Phylum	Class	Order	Family	Genus	Species
EU979477.1.1429	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Pseudorhodobacter	NA
KY770139.1.1394	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodospirillales	Thalassospiraceae	Thalassospira	NA
AB245935.1.1459	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Saprospiraceae	Aureispira	NA
KU951800.1.1271	Bacteria	Cyanobacteria	Cyanobacteriia	Synechococcales	Synechococcales_Incertae_Sedis	Schizothrix_LEGE_07164	NA
HM217045.1.1446	Bacteria	Cyanobacteria	Cyanobacteriia	Oxyphotobacteria_Incertae_Sedis	Oxyphotobacteria_Incertae_Sedis_Unknown_Family	Leptolyngbya_LEGE-06070	NA
EF654088.1.1482	Bacteria	Cyanobacteria	Cyanobacteriia	Cyanobacteriales	Geitlerinemaceae	Geitlerinema_PCC-7105	NA
KY379855.1.2514	Bacteria	Cyanobacteria	Cyanobacteriia	Cyanobacteriales	Microcystaceae	Synechocystis_PCC-6803	NA
KU951847.1.1466	Bacteria	Cyanobacteria	Cyanobacteriia	Cyanobacteriales	Cyanobacteriaceae	Symphothece_PCC-7002	NA
AB062106.1.1479	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Porphyrobacter	NA

# higher relative abundance in dark-incubated cultures than in the ocean
glm.eff.darkVSocean[glm.eff.darkVSocean$effect >= 1.5,] %>% nrow()

## [1] 18

knitr::kable(as.matrix(tax_table(biom.16s.m.genus)[glm.eff.darkVSocean[glm.eff.darkVSocean$effect >= 1.5,] %>% row.names(),]))

	Kingdom	Phylum	Class	Order	Family	Genus	Species
KC247325.1.1400	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Tropicibacter	NA
KY437089.1.1426	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Youngimonas	NA
AY546508.1.1498	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylomonadaceae_uncultured	NA
AY007296.1.1499	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylosarcina	NA
AB930603.1.1457	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylomicrobium	NA
CP002738.392605.394125	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylomonas	NA
KC963966.1.1493	Bacteria	Methanotroph	Gammaproteobacteria	Methylococcales	Methylomonadaceae	Methylobacter	NA
AB053128.1.1530	Bacteria	Proteobacteria	Gammaproteobacteria	Oceanospirillales	Alcanivoracaceae1	Alcanivorax	NA
ATTE01000001.3195237.3196755	Bacteria	Proteobacteria	Gammaproteobacteria	Thiotrichales	Thiotrichaceae	Leucothrix	NA
AB681780.1.1463	Bacteria	Proteobacteria	Gammaproteobacteria	Nitrosococcales	Methylophagaceae	Methylophaga	NA
KC534174.1.1363	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Aequorivita	NA
AB681705.1.1449	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Mesoflavibacter	NA
LN875339.1.1446	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Flavobacteriaceae	NA
HQ144198.1.1471	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Mariniflexile	NA
FJ360684.1.1470	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Meridianimaribacter	NA
KJ475509.1.1443	Bacteria	Myxococcota	Polyangia	Nannocystales	Nannocystaceae	Nannocystaceae_uncultured	NA
AB607876.1.1417	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Antarctobacter	NA
FJ232451.1.1429	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodobacterales	Rhodobacteraceae	Sedimentitalea	NA

temp.1 <- glm.eff.lightVSocean[glm.eff.lightVSocean$effect >= 1.5,] %>% rownames_to_column(var = "SpeciesID")
temp.2 <- tax_table(biom.16s.m.genus) %>% as.data.frame() %>% rownames_to_column(var = "SpeciesID")
res_sig.aldex.lightVSocean <- inner_join(temp.1, temp.2, by = "SpeciesID")
pdf(paste0(wdir, "fig/aldex-lightVSocean-scatter-genus-16s.pdf"), width = 35, height = 13)
ggplot(res_sig.aldex.lightVSocean, aes(x=reorder(Genus, diff.btw), y=diff.btw, color=Phylum)) +
    geom_point(size=3, width = 0.2) + coord_flip() +
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5)) + labs(y = "logFC", x = "Genus")
dev.off()

## png 
##   2

temp.1 <- glm.eff.darkVSocean[glm.eff.darkVSocean$effect >= 1.5, ] %>% rownames_to_column(var = "SpeciesID")
temp.2 <- tax_table(biom.16s.m.genus) %>% as.data.frame() %>% rownames_to_column(var = "SpeciesID")
res_sig.aldex.darkVSocean <- inner_join(temp.1, temp.2, by = "SpeciesID")
pdf(paste0(wdir, "fig/aldex-darkVSocean-scatter-genus-16s.pdf"), width = 35, height = 13)
ggplot(res_sig.aldex.darkVSocean, aes(x=reorder(Genus, diff.btw), y=diff.btw, color=Phylum)) +
    geom_point(size=3, width = 0.2) + coord_flip() +
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5)) + labs(y = "logFC", x = "Genus")
dev.off()

## png 
##   2

Differential abundance testing of genera with DESeq2 (minimap2)

As DESeq2 expects raw counts we will use here the minimap2 counts instead of emu counts.

#run DESeq2
#https://micca.readthedocs.io/en/latest/phyloseq.html
sample_data(biom.16s.m.genus)$Culture_setup <- as.factor(sample_data(biom.16s.m.genus)$Culture_setup)
ds <- phyloseq_to_deseq2(biom.16s.m.genus, ~ Inoculum + Culture_setup)

## converting counts to integer mode

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

ds <- DESeq(ds)

## estimating size factors
##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## final dispersion estimates

## fitting model and testing

alpha <- 0.001

# test differential expression: light.clean_gas - dark.clean_gas
res.1 <- results(ds, contrast=c("Culture_setup", "light+clean_gas", "dark+clean_gas"), alpha=alpha)
summary(res.1)

## 
## out of 1115 with nonzero total read count
## adjusted p-value < 0.001
## LFC > 0 (up)       : 35, 3.1%
## LFC < 0 (down)     : 19, 1.7%
## outliers [1]       : 15, 1.3%
## low counts [2]     : 713, 64%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

res.1 <- res.1[order(res.1$padj, na.last=NA), ]
res_sig <- res.1[(res.1$padj < alpha), ]
res_sig

## log2 fold change (MLE): Culture_setup light+clean_gas vs dark+clean_gas 
## Wald test p-value: Culture setup light.clean gas vs dark.clean gas 
## DataFrame with 54 rows and 6 columns
##                         baseMean log2FoldChange     lfcSE      stat      pvalue
##                        <numeric>      <numeric> <numeric> <numeric>   <numeric>
## HM057705.1.1470         368.2061        6.96042  0.666272  10.44683 1.51512e-25
## KC247325.1.1400          59.2219       -4.20854  0.513847  -8.19025 2.60681e-16
## KY770139.1.1394         353.2085        7.01789  0.894119   7.84895 4.19530e-15
## EF654088.1.1482        2180.9272        9.17469  1.225435   7.48688 7.05283e-14
## KJ728848.1.1458          20.0779        6.89193  0.981590   7.02119 2.19986e-12
## ...                          ...            ...       ...       ...         ...
## FJ232451.1.1429          4.59764       -2.83589  0.709342  -3.99792 6.39012e-05
## KF733610.1.1535         21.18703        4.95466  1.250006   3.96371 7.37936e-05
## AB053128.1.1530         21.91235       -3.30114  0.857355  -3.85038 1.17936e-04
## AUHJ01000047.4103.5616   3.85097        4.19135  1.095814   3.82487 1.30839e-04
## AB238790.1.1469          8.87191        3.48769  0.915463   3.80976 1.39103e-04
##                               padj
##                          <numeric>
## HM057705.1.1470        5.86352e-23
## KC247325.1.1400        5.04418e-14
## KY770139.1.1394        5.41194e-13
## EF654088.1.1482        6.82362e-12
## KJ728848.1.1458        1.70269e-10
## ...                            ...
## FJ232451.1.1429        0.000494596
## KF733610.1.1535        0.000559963
## AB053128.1.1530        0.000877715
## AUHJ01000047.4103.5616 0.000955372
## AB238790.1.1469        0.000996903

DESeq2::plotMA(res.1, alpha = alpha, ylim = c(-9,9))

res.1 %>% as.data.frame %>% filter(padj<alpha) %>% summarize(up = sum(log2FoldChange>1), down = sum(log2FoldChange<1))

##   up down
## 1 35   19

# up down
#1 35   19

res_sig <- cbind(as(res_sig, "data.frame"), as(tax_table(biom.16s.m)[rownames(res_sig), ], "matrix"))
pdf(paste0(wdir, "fig/deseq2-light_cleanVSdark_clean-scatter-genus-16s.pdf"), width = 15, height = 13)
ggplot(res_sig, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
    geom_jitter(size=3, width = 0.2) +
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
dev.off()

## png 
##   2

#View(cbind(as(res_sig, "data.frame"), as(otu_table(biom.16s.m)[rownames(res_sig), ], "matrix")))

#volcano plot
volDat <- as.data.frame(res.1)
volDat <- volDat[!is.na(volDat$padj), ]
volDat$DGE <- volDat$padj < alpha
pdf(paste0(wdir, "fig/deseq2-light_cleanVSdark_clean-volcano-genus-16s.pdf"), width = 15, height = 13)
ggplot(volDat, aes(x = log2FoldChange, y = -log10(padj), color = DGE)) + 
  geom_point(alpha = 0.50) + 
  scale_color_manual(values=c("TRUE" = "red", "FALSE" = "gray"))
dev.off

## function (which = dev.cur()) 
## {
##     if (which == 1) 
##         stop("cannot shut down device 1 (the null device)")
##     .External(C_devoff, as.integer(which))
##     dev.cur()
## }
## <bytecode: 0x5625e9ecc750>
## <environment: namespace:grDevices>

#run DESeq2
# test differential expression: light.clean_gas - light.used_gas
res.2 <- results(ds, contrast=c("Culture_setup", "light+clean_gas", "light+used_gas"), alpha=alpha)
summary(res.2)

## 
## out of 1115 with nonzero total read count
## adjusted p-value < 0.001
## LFC > 0 (up)       : 0, 0%
## LFC < 0 (down)     : 1, 0.09%
## outliers [1]       : 15, 1.3%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

res.2 <- res.2[order(res.2$padj, na.last=NA), ]
DESeq2::plotMA(res.2, alpha = alpha, ylim = c(-9,9))

res.2 %>% as.data.frame %>% filter(padj<alpha) %>% summarize(up = sum(log2FoldChange>1), down = sum(log2FoldChange<1))

##   up down
## 1  0    1

#  up down
#1  0    1
res_sig <- res.2[(res.2$padj < alpha), ]
res_sig <- cbind(as(res_sig, "data.frame"), as(tax_table(biom.16s.m.genus)[rownames(res_sig), ], "matrix"))
knitr::kable(as.matrix(res_sig %>% as.data.frame %>% filter(padj < alpha & log2FoldChange > 1) %>% dplyr::select(Genus, log2FoldChange, padj)))

Genus	log2FoldChange	padj

knitr::kable(as.matrix(res_sig %>% as.data.frame %>% filter(padj < alpha & log2FoldChange < 1) %>% dplyr::select(Genus, log2FoldChange, padj)))

	Genus	log2FoldChange	padj
JF748732.1.1472	Neptunomonas	-3.065593	0.0007264315

Differential abundance testing of genera with edgeR (minimap2)

As edgeR expects raw counts we will use here the minimap2 counts instead of emu counts.

#get genus counts
biom.16s.m.genus <- tax_glom(biom.16s.m, taxrank = "Genus")
er <- phyloseq2edgeR(biom.16s.m.genus)
#normalize
er <- edgeR::calcNormFactors(er, method = "TMM")
#metadata
meta.sub <- as.data.frame(sample_data(biom.16s.m.genus)[,c(2,3,8)])
meta.sub <- meta.sub %>% as_tibble() %>% mutate(across(.cols = 3, .fns = as.integer)) %>% mutate(across(.cols = 1:2, .fns = make.names))
#estimate dispersion
design <- model.matrix(~0 + Culture_setup + Inoculum, data=meta.sub)
er <- edgeR::estimateDisp(er, design)
plotBCV(er)

#differential expression
fit <- glmQLFit(er, design, robust = TRUE)
plotQLDisp(fit)

# test differential expression: light.clean_gas - dark.clean_gas
qlf <- glmQLFTest(fit, contrast = c(-1,1,0,0,0,0))
# create DGE table
summary(decideTests(qlf,p.value=0.001)) #FDR<0.001

##        -1*Culture_setupdark.clean_gas 1*Culture_setuplight.clean_gas
## Down                                                              11
## NotSig                                                          1089
## Up                                                                15

#       -1*Culture_setupdark.clean_gas 1*Culture_setuplight.clean_gas
#Down                                                              34
#NotSig                                                          1040
#Up                                                                41
dge.1 <- topTags(qlf, n = Inf, p.value = 0.001)$table
dge.1 %>% filter(FDR<0.001 & logFC > 1) %>% nrow()

## [1] 15

knitr::kable(as.matrix(dge.1 %>% filter(FDR<0.001 & logFC > 1) %>% dplyr::select(Genus, logFC, FDR)))

	Genus	logFC	FDR
HM057705.1.1470	Cyanobium_PCC-6307	6.170126	1.711499e-06
KY770139.1.1394	Thalassospira	6.969588	7.424660e-05
KJ728848.1.1458	Wenyingzhuangia	6.292215	8.027002e-05
CP011456.3050785.3052272	Aphanizomenon_NIES81	6.556994	8.027002e-05
EF654088.1.1482	Geitlerinema_PCC-7105	8.513339	8.027002e-05
KU951800.1.1271	Schizothrix_LEGE_07164	5.895140	9.999445e-05
JN428535.1.1439	Flexithrix	8.106133	1.046041e-04
KY379855.1.2514	Synechocystis_PCC-6803	7.440222	1.130669e-04
KY807915.1.1369	Oxyphotobacteria_Incertae_Sedis_Unknown_Family_Unknown_Family	8.295968	1.815318e-04
HM217045.1.1446	Leptolyngbya_LEGE-06070	6.937264	3.386363e-04
EU672857.1.1492	Congregibacter	4.725332	3.626087e-04
KU951847.1.1466	Symphothece_PCC-7002	7.948037	3.922224e-04
HQ674537.1.1454	Holosporaceae_uncultured	6.045458	6.340628e-04
AB062106.1.1479	Porphyrobacter	4.993965	6.700702e-04
GQ441334.1.1453	Spirosomaceae_uncultured	4.724273	7.052793e-04

dge.1 %>% filter(FDR<0.001 & logFC < 1) %>% nrow()

## [1] 11

knitr::kable(as.matrix(dge.1 %>% filter(FDR<0.001 & logFC < 1) %>% dplyr::select(Genus, logFC, FDR)))

	Genus	logFC	FDR
CP018572.2205903.2207370	Marivivens	-7.271044	2.249065e-05
KC247325.1.1400	Tropicibacter	-4.721854	2.249065e-05
EF471465.1.1394	Litorimicrobium	-6.666238	8.027002e-05
KC963966.1.1493	Methylobacter	-7.669797	8.027002e-05
FJ869046.1.1459	Aestuariicoccus	-6.639487	8.027002e-05
ATTE01000001.3195237.3196755	Leucothrix	-5.217210	8.027002e-05
DQ234217.1.1519	HIMB11	-6.301259	1.167209e-04
CP022207.429135.430671	Francisella	-7.170323	1.656684e-04
AB550558.1.1436	Primorskyibacter	-4.244944	3.050166e-04
AB681705.1.1449	Mesoflavibacter	-5.602549	7.052793e-04
JX844498.1.1439	Subsaxibacter	-4.737066	7.052793e-04

res_sig <- dge.1 %>% filter(FDR<0.001)
pdf(paste0(wdir, "fig/edger-light_cleanVSdark_clean-scatter-genus-16s.pdf"), width = 15, height = 13)
ggplot(res_sig, aes(x=Genus, y=logFC, color=Phylum)) +
    geom_jitter(size=3, width = 0.2) +
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
dev.off()

## png 
##   2

# test differential expression: (light.clean_gas + light.used_gas) - dark.clean_gas
qlf <- glmQLFTest(fit, contrast = c(-1,1/2,1/2,0,0,0))
# create DGE table
summary(decideTests(qlf,p.value=0.001)) #FDR<0.001

##        -1*Culture_setupdark.clean_gas 0.5*Culture_setuplight.clean_gas 0.5*Culture_setuplight.used_gas
## Down                                                                                                38
## NotSig                                                                                            1053
## Up                                                                                                  24

#       -1*Culture_setupdark.clean_gas 0.5*Culture_setuplight.clean_gas 0.5*Culture_setuplight.used_gas
#Down                                                                                                64
#NotSig                                                                                             997
#Up                                                                                                  54
dge.2 <- topTags(qlf, n = Inf, p.value = 0.001)$table
dge.2 %>% filter(FDR<0.001 & logFC > 1) %>% nrow()

## [1] 24

knitr::kable(as.matrix(dge.2 %>% filter(FDR<0.001 & logFC > 1) %>% dplyr::select(Genus, logFC, FDR)))

	Genus	logFC	FDR
HM057705.1.1470	Cyanobium_PCC-6307	5.258539	2.949909e-06
KY770139.1.1394	Thalassospira	7.064298	1.402626e-05
KJ728848.1.1458	Wenyingzhuangia	6.417309	1.445156e-05
CP011456.3050785.3052272	Aphanizomenon_NIES81	6.595324	1.838711e-05
KU951800.1.1271	Schizothrix_LEGE_07164	5.881805	2.783457e-05
EF654088.1.1482	Geitlerinema_PCC-7105	8.626385	2.877699e-05
KY807915.1.1369	Oxyphotobacteria_Incertae_Sedis_Unknown_Family_Unknown_Family	8.380365	6.811987e-05
HM217045.1.1446	Leptolyngbya_LEGE-06070	7.071307	1.166464e-04
KU951847.1.1466	Symphothece_PCC-7002	8.111378	1.429601e-04
KU578816.1.1368	NS2b_marine_group	4.970993	1.590979e-04
GU230385.1.1498	Peredibacter	4.566924	2.417101e-04
AB062106.1.1479	Porphyrobacter	4.908613	2.492908e-04
JN428535.1.1439	Flexithrix	6.724950	2.614798e-04
KY379855.1.2514	Synechocystis_PCC-6803	6.255655	2.903490e-04
LC275989.1.1406	Limnothrix	7.158101	3.835265e-04
EU672857.1.1492	Congregibacter	4.295788	3.903955e-04
DQ015830.1.1494	Oleiphilus	3.914621	3.928785e-04
KM587636.1.1505	Cyclobacteriaceae	4.888703	6.036249e-04
KM019986.1.1462	Pleurocapsa_PCC-7319	7.026484	6.969247e-04
ALVO01000007.393816.395283	Geminocystis_PCC-6308	5.807228	7.004955e-04
KF733610.1.1535	Gimesia	4.068907	8.445897e-04
KM820575.1.1347	Paraglaciecola	3.324822	8.751609e-04
AB245935.1.1459	Aureispira	7.677378	8.751609e-04
AM279196.8.1489	NS11-12_marine_group	3.984948	8.845655e-04

dge.2 %>% filter(FDR<0.001 & logFC < 1) %>% nrow()

## [1] 38

knitr::kable(as.matrix(dge.2 %>% filter(FDR<0.001 & logFC < 1) %>% dplyr::select(Genus, logFC, FDR)))

	Genus	logFC	FDR
KC963966.1.1493	Methylobacter	-7.662905	2.635192e-07
KC247325.1.1400	Tropicibacter	-4.705569	2.635192e-07
CP022207.429135.430671	Francisella	-8.211281	1.058575e-06
EF471465.1.1394	Litorimicrobium	-6.611537	1.434257e-06
ATTE01000001.3195237.3196755	Leucothrix	-5.029577	1.984477e-06
CP018572.2205903.2207370	Marivivens	-6.066694	1.984477e-06
AB681705.1.1449	Mesoflavibacter	-5.956536	5.265970e-06
AB930603.1.1457	Methylomicrobium	-7.344289	1.092613e-05
KC534174.1.1363	Aequorivita	-5.631045	1.387040e-05
FJ869046.1.1459	Aestuariicoccus	-5.424786	1.445156e-05
AB500095.1.1519	Marinicella	-4.125019	1.445156e-05
AY546508.1.1498	Methylomonadaceae_uncultured	-7.765463	1.445156e-05
KY437089.1.1426	Youngimonas	-5.755412	1.445156e-05
JX844498.1.1439	Subsaxibacter	-4.858413	1.668599e-05
FJ264570.1.1500	pItb-vmat-59	-6.536707	4.851993e-05
DQ234217.1.1519	HIMB11	-4.806107	4.851993e-05
FMJP01000825.542.2072	Methyloprofundus	-5.165890	4.851993e-05
AY007296.1.1499	Methylosarcina	-7.028768	4.851993e-05
CP015231.476529.477997	Ruegeria	-3.581767	4.851993e-05
AB550558.1.1436	Primorskyibacter	-3.659586	6.277639e-05
KF465276.1.1344	Methylomarinum	-5.917652	9.508996e-05
AB453959.1.1457	Marine_Methylotrophic_Group_2	-6.184360	9.692300e-05
KC290421.1.1525	Crenothrix	-6.180702	1.187234e-04
LN875339.1.1446	Flavobacteriaceae	-3.820729	1.252776e-04
AY661691.1.1473	Tenacibaculum	-3.839027	1.252776e-04
KF437680.1.1446	Vitellibacter	-4.426200	1.320965e-04
FJ712568.1.1394	Milano-WF1B-03	-5.738365	1.590979e-04
AB607876.1.1417	Antarctobacter	-3.455238	2.375354e-04
AB053128.1.1530	Alcanivorax	-3.793586	2.375354e-04
MKSF01000014.384003.385521	Taibaiella	-5.627458	2.479392e-04
CP019307.2493152.2494620	Phaeobacter	-3.596651	2.479392e-04
GU451533.1.1243	Schleiferia	-4.788819	2.479392e-04
EU287117.1.1428	Albimonas	-3.411079	4.003033e-04
JQ197367.1.1350	Aquibacter	-2.789719	4.676467e-04
HQ144198.1.1471	Mariniflexile	-3.665439	5.352124e-04
FJ436725.1.1370	Celeribacter	-2.747662	7.774477e-04
JX412960.1.1486	Flavirhabdus	-4.582240	8.335923e-04
JX391335.1.1502	Methylophagaceae_uncultured	-4.924455	9.761934e-04

pdf(paste0(wdir, "fig/edger-light_clean_usedVSdark_clean-scatter-genus-16s.pdf"), width = 15, height = 13)
ggplot(dge.2, aes(x=Genus, y=logFC, color=Phylum)) +
    geom_jitter(size=3, width = 0.2) +
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
dev.off()

## png 
##   2

# test differential expression: light.clean_gas - light.used_gas
qlf <- glmQLFTest(fit, contrast = c(0,1,-1,0,0,0))
# create DGE table
summary(decideTests(qlf,p.value=0.001)) #FDR<0.001

##        1*Culture_setuplight.clean_gas -1*Culture_setuplight.used_gas
## Down                                                               1
## NotSig                                                          1114
## Up                                                                 0

#       1*Culture_setuplight.clean_gas -1*Culture_setuplight.used_gas
#Down                                                               4
#NotSig                                                          1098
#Up                                                                 0
dge.3 <- topTags(qlf, n = Inf)$table
dge.3 %>% filter(FDR<0.001 & logFC > 1) %>% nrow()

## [1] 0

knitr::kable(as.matrix(dge.3 %>% filter(FDR<0.001 & logFC > 1) %>% dplyr::select(Genus, logFC, FDR)))

Genus	logFC	FDR

dge.3 %>% filter(FDR<0.001 & logFC < 1) %>% nrow()

## [1] 1

knitr::kable(as.matrix(dge.3 %>% filter(FDR<0.001 & logFC < 1) %>% dplyr::select(Genus, logFC, FDR)))

	Genus	logFC	FDR
AM279196.8.1489	NS11-12_marine_group	-5.153853	0.0009801917

res_sig <- dge.3 %>% filter(FDR<0.001)
pdf(paste0(wdir, "fig/edger-light_cleanVSlight_used-scatter-genus.pdf"), width = 15, height = 13)
ggplot(res_sig, aes(x=Genus, y=logFC, color=Phylum)) +
    geom_jitter(size=3, width = 0.2) +
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
dev.off()

## png 
##   2

# test differential expression: ocean - (light.clean_gas + light.used_gas + dark.clean_gas)/3
qlf <- glmQLFTest(fit, contrast = c(-1/3,-1/3,-1/3,1,0,0))
# create DGE table
summary(decideTests(qlf,p.value=0.001)) #FDR<0.001

##        -0.333333333333333*Culture_setupdark.clean_gas -0.333333333333333*Culture_setuplight.clean_gas -0.333333333333333*Culture_setuplight.used_gas 1*Culture_setupocean
## Down                                                                                                                                                                   23
## NotSig                                                                                                                                                                934
## Up                                                                                                                                                                    158

#       -0.333333333333333*Culture_setupdark.clean_gas -0.333333333333333*Culture_setuplight.clean_gas -0.333333333333333*Culture_setuplight.used_gas 1*Culture_setupocean
#Down                                                               45
#NotSig                                                            828
#Up                                                                242
dge.4 <- topTags(qlf, n = Inf, p.value = 0.001)$table
dge.4 %>% filter(FDR<0.001 & logFC > 1) %>% nrow()

## [1] 158

knitr::kable(as.matrix(dge.4 %>% filter(FDR<0.001 & logFC > 1) %>% dplyr::select(Genus, logFC, FDR)))

	Genus	logFC	FDR
AY548995.1.1382	Lentimicrobiaceae	4.299767	3.128146e-10
EU799083.1.1494	NS5_marine_group	4.588771	6.742205e-10
KF741413.1.1479	Bacteroidetes_VC2.1_Bac22	4.702386	6.742205e-10
KT266595.1.1504	Desulforhopalus	4.164808	1.796639e-09
KC425506.1.1449	SAR116_clade	4.700587	1.796639e-09
DQ234096.1.1547	Arcobacteraceae_uncultured	4.398214	2.280214e-09
CP003984.2903910.2905362	Planktomarina	4.363532	2.280214e-09
JX391370.1.1507	Thiomicrorhabdus	4.352153	2.280214e-09
JN977181.1.1499	Sulfurovum	4.311595	4.011571e-09
DQ351806.1.1480	Bacteroidetes_BD2-2	4.161159	4.844382e-09
EU432455.1.1361	Marinifilaceae_uncultured	4.431837	4.890420e-09
KX172972.1.1444	Draconibacterium	3.926086	1.741132e-08
AM279194.8.1485	Crocinitomicaceae_uncultured	3.959948	1.741132e-08
CP000084.511358.512831	Clade_Ia	4.471599	1.741132e-08
DQ234254.1.1547	Pseudarcobacter	4.296199	1.997420e-08
KX177096.1.1491	Lachnospiraceae_uncultured	4.064750	1.997420e-08
JN107026.1.1478	Fermentibacteraceae	5.468135	1.997420e-08
DQ677856.1.1392	Prolixibacteraceae	4.428966	1.997420e-08
KT266596.1.1481	Marinifilum	4.209685	2.499360e-08
FJ545595.1.1521	NS4_marine_group	4.007211	4.274081e-08
JX391804.1.1475	Sulfurimonas	3.964505	4.274081e-08
AACY020462030.661.2167	Candidatus_Actinomarina	4.015651	4.274081e-08
HM057617.1.1504	SAR86_clade	3.997552	4.274081e-08
FR733673.1.1549	Desulfobacter	3.960248	4.274081e-08
KF268736.1.1370	Roseimarinus	4.784350	4.355030e-08
KC492879.1.1435	MSBL3	4.466855	5.493685e-08
EF632644.1.1398	Hungateiclostridiaceae_uncultured	4.693646	5.904163e-08
KF650753.1.1478	Halarcobacter	4.147976	6.730820e-08
FR684261.1.1500	Pseudohongiella	4.120351	6.879621e-08
JF344182.1.1488	Marinilabiliaceae_uncultured	3.733532	7.015661e-08
KR825117.1.1534	R76-B128	3.964217	8.377625e-08
FJ672721.1.1394	Alkaliflexus	4.446077	1.153754e-07
AY177803.1.1522	Desulfoconvexum	4.175328	1.153754e-07
GQ274297.1.1512	Ichthyenterobacterium	3.333161	1.153754e-07
HQ008573.1.1509	hgcI_clade	4.161558	1.198956e-07
KX097755.1.1517	MSBL8	3.982331	1.359139e-07
HM234096.1.1443	Lutibacter	3.703584	1.654225e-07
FJ202580.1.1470	Defluviitaleaceae_UCG-011	4.014297	1.849523e-07
JQ580521.1.1521	Desulfosarcinaceae_uncultured	4.012735	2.083910e-07
JF344429.1.1502	Spirochaeta_2	4.044421	2.113200e-07
AACY020059431.1419.2884	Clade_III	3.942214	2.116204e-07
KC582602.1.1415	Fluviicola	3.105267	2.129226e-07
FLOH01001865.2173.3719	Desulfatiglans	3.876941	2.582172e-07
JF344214.1.1499	Sphaerochaeta	3.948439	3.244416e-07
FJ436732.1.1375	Lentibacter	4.034904	3.597343e-07
HM057738.1.1543	OM43_clade	4.014748	3.826552e-07
GQ356973.1.1496	livecontrolB21	3.813504	4.843966e-07
GQ246314.1.1499	Woeseia	4.527468	5.946380e-07
KF799284.1.1451	Synechococcus_CC9902	3.298238	6.611743e-07
FM204974.1.1427	Dysgonomonadaceae_uncultured	4.257319	8.567944e-07
JX391453.1.1506	Desulfocapsaceae_uncultured	4.039990	9.585658e-07
AACY020291711.299.1802	Marinobacterium	3.882655	1.051261e-06
KX088587.1.1518	Sediminispirochaeta	4.178136	1.765495e-06
JX016104.1.1480	Flavicella	5.492827	2.218415e-06
HQ203800.1.1480	Arcobacteraceae	4.094904	2.226069e-06
AACY020468607.261.1725	Clade_II	3.872826	2.630271e-06
EU358688.1.1512	Treponema	3.751128	3.550629e-06
FR684473.1.1498	RS62_marine_group	4.057283	3.633471e-06
EU795171.11703.13180	Thiotrichaceae_uncultured	3.665790	4.549235e-06
CP017259.618357.619876	Formosa	3.161721	4.649038e-06
EU239499.1.1482	Lutimonas	3.031947	6.045904e-06
EU795239.16117.17616	Ulvibacter	2.908831	6.983980e-06
JX403044.1.1509	Arcobacter	4.568130	6.983980e-06
AB546078.1.1345	Anaerovorax	4.685314	6.983980e-06
JX391399.1.1481	Saccharicrinis	3.974172	9.504767e-06
EF061973.1.1483	Labilibacter	3.896357	1.260831e-05
FPLK01002149.16.1518	CL500-29_marine_group	3.969527	1.421805e-05
AB630524.1.1448	Bacteroidetes_vadinHA17	4.053649	1.819234e-05
JQ925053.1.1549	Candidatus_Omnitrophus	4.745160	1.819234e-05
GQ249615.1.1406	Latescibacterota	3.734224	1.973933e-05
EF044234.1.1421	Roseibacterium	2.309185	2.074793e-05
AY771944.1.1521	Desulfosarcina	4.185723	2.074793e-05
EU864446.1.1617	[Eubacterium]_coprostanoligenes_group	3.813844	2.198103e-05
HQ405605.1.1522	SG8-4	4.557494	2.209915e-05
EU642572.1.1399	Thiothrix	3.593740	2.224185e-05
JX225893.1.1480	WCHB1-32	4.058192	2.271422e-05
DQ811934.1.1510	Calditrichaceae	4.947964	3.020438e-05
HE652874.1.1297	Lentimicrobium	4.029582	3.036462e-05
KX097528.1.1498	Caldithrix	3.981808	3.156894e-05
GQ356959.1.1520	Sva0485	4.461833	3.263169e-05
EU287304.1.1489	Eudoraea	3.744319	3.360131e-05
AY344418.1.1493	NS9_marine_group	2.779503	3.369507e-05
CP019352.2140807.2142329	Lacinutrix	3.982231	3.650767e-05
LC192989.1.1290	Stappiaceae	4.872698	4.300884e-05
KR824986.1.1515	Aminicenantales	3.586835	4.410536e-05
HM057716.1.1496	SAR324_clade(Marine_group_B)	3.791822	5.296689e-05
AB630756.1.1458	Arenicellaceae_uncultured	4.418655	6.308305e-05
CP021023.969900.971425	L21-RPul-D3	4.613935	6.308305e-05
FJ517016.1.1505	MidBa8	4.946262	7.698632e-05
DQ831540.1.1516	Desulfobacteraceae_uncultured	4.063789	7.698632e-05
FR685242.1.1486	NS3a_marine_group	3.474037	7.698632e-05
JX017163.1.1500	MB11C04_marine_group	4.364665	9.022946e-05
FM242394.1.1450	Planktothricoides_SR001	4.262996	9.022946e-05
HG971004.1.1492	Rikenellaceae_RC9_gut_group	3.474231	9.022946e-05
FM242289.1.1473	Leptotrichiaceae_uncultured	3.631345	9.192667e-05
EU801277.1.1280	SAR11_clade_uncultured_uncultured	3.853660	1.116840e-04
JF514267.1.1484	Christensenellaceae_uncultured	3.639531	1.131463e-04
EU592362.1.1510	Marinimicrobia_(SAR406_clade)	4.201891	1.131463e-04
DQ811912.1.1485	Salinirepens	4.054921	1.141285e-04
JN119244.1.1494	MWH-UniP1_aquatic_group	4.133162	1.141368e-04
KX933250.1.1378	PeM15	2.796338	1.176452e-04
HM231145.1.1450	Ruminococcus	4.843789	1.347174e-04
AY188291.1.1480	Christensenellaceae_R-7_group	4.396753	1.354872e-04
JQ579954.1.1492	IheB3-7	3.343783	1.365589e-04
AY337519.1.1512	Clostridium_sensu_stricto_1	3.275686	1.608877e-04
KR825137.1.1546	JS1	3.954690	1.619184e-04
HQ405627.1.1487	Demequina	4.241660	1.619184e-04
KJ752776.1.1445	Vicingus	2.924077	1.619184e-04
EU719657.1.1388	Dethiosulfovibrio	4.870964	1.619184e-04
AJ431235.1.1477	Lenti-02	3.801318	1.695557e-04
FR683698.1.1427	Octadecabacter	3.049596	1.807085e-04
EU488104.1.1522	BD2-11_terrestrial_group	4.624778	1.810754e-04
FO203503.30555.32102	Desulfobacula	4.229121	1.882647e-04
AB176674.1.1504	Lishizhenia	4.466516	2.003596e-04
FJ825870.1.1488	Cryomorphaceae_uncultured	3.277393	2.020620e-04
LC193132.1.1465	Cocleimonas	3.395675	2.067887e-04
JX391227.1.1480	MAT-CR-H4-C10	4.464705	2.074201e-04
GQ246449.1.1491	4572-13	3.466848	2.105930e-04
HQ132402.1.1485	LCP-89	3.421611	2.367440e-04
HQ326320.1.1506	Rhodanobacter	4.234875	2.392887e-04
KP091115.1.1623	SEEP-SRB1	4.388404	2.447958e-04
KX097751.1.1498	UCG-012	4.593284	2.465312e-04
JQ197367.1.1350	Aquibacter	2.350957	2.555607e-04
HM992535.1.1347	Thermovirga	3.686608	2.555607e-04
AY177796.1.1501	[Desulfobacterium]_catecholicum_group	4.579529	2.818517e-04
FJ664782.1.1404	Paludibacteraceae_uncultured	4.387308	2.818517e-04
JN458997.1.1446	0319-6G20	4.601291	2.877755e-04
KC665950.1.1514	Anaerophaga	4.072531	2.919247e-04
HF922341.1.1384	Dethiosulfatibacter	4.369118	3.109274e-04
JN637820.1.1471	Sulfurospirillum	4.307876	3.236673e-04
GQ383925.1.1483	Algoriella	3.992145	3.313441e-04
GU196238.1.1492	Paludibacter	3.913049	3.313441e-04
JF268397.1.1477	SB-5	4.238587	3.611305e-04
AY548942.1.1471	Anaerolineaceae_uncultured	4.295069	3.715710e-04
LS997911.1.1416	Desulfovibrio	4.560066	3.852181e-04
DQ302106.1.1491	Coraliomargarita	4.388408	4.172618e-04
FN433319.1.1575	Leeuwenhoekiella	4.045886	4.247481e-04
JX391727.1.1479	Prolixibacteraceae_uncultured	3.927565	5.049171e-04
KM255690.1.1435	Thioclava	2.756104	5.092086e-04
KR825155.1.1455	Candidatus_Thiobios	3.110537	5.092086e-04
EU799705.1.1457	PS1_clade	3.029088	5.452319e-04
JX391178.1.1531	Desulfitobacteriales_uncultured_uncultured	3.813120	5.467212e-04
EF574272.1.1491	PB19	2.464247	5.615888e-04
KF741432.1.1496	Clostridia_vadinBB60_group	3.959153	5.790300e-04
AB212893.1.1508	Luteibacter	4.284996	5.907625e-04
EF400726.1.1488	Odoribacter	4.200624	5.940350e-04
EF655677.1.1491	Guggenheimella	3.603800	7.171432e-04
KF799120.1.1498	HOC36	3.739430	7.447251e-04
FJ189532.1.1479	SBR1031	3.979767	7.447251e-04
HQ891651.1.1387	Planktotalea	2.551894	7.817720e-04
JX016885.1.1476	Flavobacteriaceae_uncultured	2.164628	8.136996e-04
GQ850575.1.1521	Subgroup_23	3.624799	8.136996e-04
JN496611.1.1463	LD1-PA32	4.082652	9.179859e-04
JN387332.1.1319	RBG-13-54-9	3.455901	9.179859e-04
GU451382.1.1229	Rubidimonas	3.267613	9.179859e-04
JQ580392.1.1467	Acetobacterium	4.135390	9.707469e-04
FJ960414.1.1392	Weeksellaceae_uncultured	3.665669	9.836180e-04
EU802139.1.1502	OM182_clade	3.419255	9.836180e-04

dge.4 %>% filter(FDR<0.001 & logFC < 1) %>% nrow()

## [1] 23

knitr::kable(as.matrix(dge.4 %>% filter(FDR<0.001 & logFC < 1) %>% dplyr::select(Genus, logFC, FDR)))

	Genus	logFC	FDR
MH197117.1.1432	Nioella	-3.534586	9.504767e-06
KJ754648.1.1322	Nautella	-5.488589	9.504767e-06
JX426065.1.1385	Taeseokella	-4.372977	9.785866e-06
MH396689.1.1476	Cohaesibacter	-4.197273	1.156314e-05
AB627076.1.1425	Lutimaribacter	-4.900967	2.007194e-05
MK053889.1.1442	Marinovum	-4.118054	2.195942e-05
CP002738.392605.394125	Methylomonas	-5.006456	2.567470e-05
GU324080.1.1459	Marivita	-5.117043	2.677103e-05
EU221274.1.1422	Seohaeicola	-3.409820	3.199306e-05
AB062106.1.1479	Porphyrobacter	-4.333170	7.904373e-05
CP012040.3233347.3234858	Cyclobacterium	-4.944991	9.106964e-05
AB238790.1.1469	Aliiglaciecola	-4.188013	1.141368e-04
AJ786600.1.1449	Hoeflea	-3.567957	1.204804e-04
EU979477.1.1429	Pseudorhodobacter	-3.648229	1.365589e-04
KC160833.1.1501	Marinobacter	-3.195230	1.996241e-04
FJ624296.1.1434	Sagittula	-2.901668	2.067887e-04
KC120671.1.1353	Sneathiella	-3.699667	2.288779e-04
GQ441325.1.1482	Lewinella	-4.127695	2.887393e-04
HQ316943.1.1422	Oricola	-3.100572	4.714132e-04
JN119141.1.1383	Rhodobacteraceae	-2.515310	5.574295e-04
KP137426.1.1396	Denitromonas	-3.585768	7.447251e-04
MH463965.1.1530	Mesorhizobium	-3.616027	7.845981e-04
JX223093.1.1487	Leucobacter	-3.626221	8.136996e-04

res_sig <- dge.4 %>% filter(FDR<0.001)
pdf(paste0(wdir, "fig/edger-oceanVSinoculum-scatter-genus-16s.pdf"), width = 35, height = 13)
ggplot(res_sig, aes(x=Genus, y=logFC, color=Phylum)) +
    geom_jitter(size=3, width = 0.2) +
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
dev.off()

## png 
##   2

Linear time course trends under light of genera with edgeR (minimap2)

# time course trend analysis
library(splines)
meta.light <- meta.sub %>% filter(Culture_setup == "light.clean_gas" | Culture_setup == "light.used_gas" | Culture_setup == "ocean") %>% filter(Inoculum != "Hornbaek.shallow")
X <- ns(meta.light$Time_point, df=1)
design.light <- model.matrix(~X + Inoculum, data=meta.light)
biom.16s.m.light <- prune_samples(sample_names(biom.16s.m.genus) %in% c("flaregas-1","flaregas-2","flaregas-3","flaregas-10","flaregas-11","flaregas-12","flaregas-13","flaregas-14","flaregas-15","flaregas-16","flaregas-17","flaregas-18","flaregas-19","flaregas-20"), biom.16s.m.genus)
er <- phyloseq2edgeR(biom.16s.m.light)
er <- edgeR::calcNormFactors(er, method = "TMM")
er <- edgeR::estimateDisp(er, design.light)
plotBCV(er)

sqrt(er$common.dispersion)

## [1] 0.9915841

fit <- glmQLFit(er, design.light, robust = TRUE)
plotQLDisp(fit)

qlf <- glmQLFTest(fit, coef=2)
summary(decideTests(qlf, p.value=0.01))

##           X
## Down     85
## NotSig 1002
## Up       28

#          X
#Down     60
#NotSig 1033
#Up       22

#filter genera that where detected with Emu (filter low abundant genera)
dge.5 <- topTags(qlf, n = Inf, p.value = 1)$table
dge.5 <- dge.5[biom.16s.genus.minimap2emu$SpeciesID,]

# visualize the fitted spline curves for the top four enriched species under cultivation conditions
dge.5.u <- dge.5 %>% filter(FDR<0.01 & logFC > 1)
dge.5.u %>% nrow()

## [1] 21

logCPM.obs <- cpm(er, log=TRUE, prior.count=qlf$prior.count)
logCPM.fit <- cpm(qlf, log=TRUE)
pdf(paste0(wdir, "fig/edger-time-course-light-genus-16s-1.pdf"), width = 15, height = 15)
par(mfrow=c(3,7))
for(i in 1:nrow(dge.5.u)) {
 FlybaseID <- row.names(dge.5.u)[i]
 Symbol <- dge.5.u$Genus[i]
 logCPM.obs.i <- logCPM.obs[FlybaseID,]
 logCPM.fit.i <- logCPM.fit[FlybaseID,]
 plot(meta.light$Time_point, logCPM.obs.i, ylab="log-CPM", xlab="Time point", main=Symbol, pch=16)
 lines(meta.light$Time_point, logCPM.fit.i, col="red", lwd=2)
}
dev.off()

## png 
##   2

Linear time course trends in the dark of genera with edgeR (minimap2)

# time course trend analysis
library(splines)
meta.dark <- meta.sub %>% filter(Culture_setup == "dark.clean_gas" | Culture_setup == "ocean") %>% filter(Inoculum != "Helsingor")
X <- ns(meta.dark$Time_point, df=1)
design.dark <- model.matrix(~X + Inoculum, data=meta.dark)
biom.16s.m.dark <- prune_samples(sample_names(biom.16s.m.genus) %in% c("flaregas-11","flaregas-12","flaregas-21","flaregas-22","flaregas-23","flaregas-24","flaregas-25","flaregas-26","flaregas-27","flaregas-28","flaregas-29","flaregas-30"), biom.16s.m.genus)
er <- phyloseq2edgeR(biom.16s.m.dark)
er <- edgeR::calcNormFactors(er, method = "TMM")
er <- edgeR::estimateDisp(er, design.dark)
plotBCV(er)

sqrt(er$common.dispersion)

## [1] 0.7714087

fit <- glmQLFit(er, design.dark, robust = TRUE)
plotQLDisp(fit)

qlf <- glmQLFTest(fit, coef=2)
summary(decideTests(qlf, p.value=0.01))

##          X
## Down   233
## NotSig 833
## Up      49

#         X
#Down   248
#NotSig 837
#Up      30

#filter genera that where detected with Emu (filter low abundant genera)
dge.6 <- topTags(qlf, n = Inf, p.value = 1)$table
dge.6 <- dge.6[biom.16s.genus.minimap2emu$SpeciesID,]

# visualize the fitted spline curves for the top four enriched species under cultivation conditions
dge.6.u <- dge.6 %>% filter(FDR<0.01 & logFC > 1)
dge.6.u %>% nrow()

## [1] 47

logCPM.obs <- cpm(er, log=TRUE, prior.count=qlf$prior.count)
logCPM.fit <- cpm(qlf, log=TRUE)
pdf(paste0(wdir, "fig/edger-time-course-dark-genus-16s-1.pdf"), width = 15, height = 15)
par(mfrow=c(4,7))
for(i in 1:28) {
 FlybaseID <- row.names(dge.6.u)[i]
 Symbol <- dge.6.u$Genus[i]
 logCPM.obs.i <- logCPM.obs[FlybaseID,]
 logCPM.fit.i <- logCPM.fit[FlybaseID,]
 plot(meta.dark$Time_point, logCPM.obs.i, ylab="log-CPM", xlab="Time point", main=Symbol, pch=16)
 lines(meta.dark$Time_point, logCPM.fit.i, col="red", lwd=2)
}
dev.off()

## png 
##   2

pdf(paste0(wdir, "fig/edger-time-course-dark-genus-16s-2.pdf"), width = 15, height = 15)
par(mfrow=c(3,7))
for(i in 29:nrow(dge.6.u)) {
 FlybaseID <- row.names(dge.6.u)[i]
 Symbol <- dge.6.u$Genus[i]
 logCPM.obs.i <- logCPM.obs[FlybaseID,]
 logCPM.fit.i <- logCPM.fit[FlybaseID,]
 plot(meta.dark$Time_point, logCPM.obs.i, ylab="log-CPM", xlab="Time point", main=Symbol, pch=16)
 lines(meta.dark$Time_point, logCPM.fit.i, col="red", lwd=2)
}
dev.off()

## png 
##   2

Compare results of DESeq2 and edgeR on genus level

alpha <- 0.001
# test differential expression: light.clean_gas - dark.clean_gas
gene_list <- list(
  DESeq2 = res.1 %>% as.data.frame %>% dplyr::filter(padj < alpha & log2FoldChange<1) %>% rownames,
  edgeR = dge.1 %>% dplyr::filter(FDR < alpha & logFC<1) %>% rownames
)
#venn.diagram(x = gene_list, filename = "venn.light-cleanVSdark-clean.down.tiff")
p.a <- venn(gene_list)

gene_list <- list(
  DESeq2 = res.1 %>% as.data.frame %>% dplyr::filter(padj < alpha & log2FoldChange>1) %>% rownames,
  edgeR = dge.1 %>% dplyr::filter(FDR < alpha & logFC>1) %>% rownames
)
p.b <- venn(gene_list)

Compare light_clean_usedVSdark results of ALDEX and edgeR on genus level

alpha <- 0.001
# test differential expression: (light.clean_gas + light.used_gas) - dark.clean_gas
gene_list <- list(
  aldex = glm.eff.lightVSdark[glm.eff.lightVSdark$effect >= 1,] %>% row.names(), 
  edgeR = dge.2 %>% filter(FDR < alpha & logFC > 1) %>% rownames
)
p.a <- venn(gene_list)

knitr::kable(as.matrix(tax_table(biom.16s.m.genus)[attr(p.a,"intersections")$aldex,]))
gene_list <- list(
  aldex = glm.eff.lightVSdark[glm.eff.lightVSdark$effect <= -1,] %>% row.names(), 
  edgeR = dge.2 %>% filter(FDR < alpha & logFC < 1) %>% rownames
)
p.a <- venn(gene_list)

knitr::kable(as.matrix(tax_table(biom.16s.m.genus)[attr(p.a,"intersections")$aldex,]))

Compare results of ocean VS inoculum and time course trends in edgeR on genus level

alpha <- 0.01
# test differential expression: ocean - (light.clean_gas + light.used_gas + dark.clean_gas)/3
#light
gene_list <- list(
  edgeR.contrast = dge.4 %>% filter(FDR < alpha & logFC > 1) %>% rownames,
  edgeR.splines = dge.5 %>% filter(FDR < alpha & logFC < 1) %>% rownames
)
p.a <- venn(gene_list)

gene_list <- list(
  edgeR.contrast = dge.4 %>% filter(FDR < alpha & logFC < 1) %>% rownames,
  edgeR.splines = dge.5 %>% filter(FDR < alpha & logFC > 1) %>% rownames
)
p.b <- venn(gene_list)

#dark
gene_list <- list(
  edgeR.contrast = dge.4 %>% filter(FDR < alpha & logFC > 1) %>% rownames,
  edgeR.splines = dge.6 %>% filter(FDR < alpha & logFC < 1) %>% rownames
)
p.a <- venn(gene_list)

gene_list <- list(
  edgeR.contrast = dge.4 %>% filter(FDR < alpha & logFC < 1) %>% rownames,
  edgeR.splines = dge.6 %>% filter(FDR < alpha & logFC > 1) %>% rownames
)
p.b <- venn(gene_list)

Network construction and characterization

Semi-Parametric Rank-based approach for INference in Graphical model (SPRING)

#https://github.com/stefpeschel/NetCoMi
library(NetCoMi)

## Loading required package: SpiecEasi

## 
## Attaching package: 'SpiecEasi'

## The following object is masked from 'package:MASS':
## 
##     fitdistr

## 

net_spring <- netConstruct(biom.16s,
                           filtTax = "highestFreq",
                           filtTaxPar = list(highestFreq = 50),
                           measure = "spring",
                           measurePar = list(nlambda=50, 
                                             rep.num=50,
                                             Rmethod = "approx"),
                           normMethod = "none", 
                           zeroMethod = "none",
                           sparsMethod = "none", 
                           dissFunc = "signed",
                           verbose = 3,
                           seed = 123456)

## Checking input arguments ...

## Done.

## Data filtering ...

## 928 taxa removed.

## 50 taxa and 24 samples remaining.

## 
## Calculate 'spring' associations ... 
## Registered S3 method overwritten by 'dendextend':
##   method     from 
##   rev.hclust vegan
## Registered S3 method overwritten by 'seriation':
##   method         from 
##   reorder.hclust vegan
## Done.

#Adjacency values computed using the "signed" transformation (values different from 0 and 1 will be edges in the network):
hist(net_spring$assoMat1, 100, xlim = c(-1, 1), ylim = c(0, 400),
     xlab = "Estimated correlation", 
     main = "Estimated correlations after sparsification")

net_spring_unsigned <- netConstruct(data = net_spring$assoEst1,
                           dataType = "correlation",
                           sparsMethod = "threshold",
                           thresh = 0.3,
                           dissFunc = "unsigned",
                           verbose = 3)

## Checking input arguments ... Done.
## 
## Sparsify associations via 'threshold' ... Done.

#Adjacency values computed using the "unsigned" transformation:
hist(net_spring_unsigned$adjaMat1, 100, ylim = c(0, 400),
     xlab = "Adjacency values", 
     main = "Adjacencies (with \"unsigned\" transformation)")

temp.nw <- net_spring
props_spring <- netAnalyze(temp.nw, 
                           centrLCC = TRUE,
                           clustMethod = "cluster_fast_greedy",
                           hubPar = "eigenvector",
                           weightDeg = FALSE, normDeg = FALSE)

#?summary.microNetProps
print(summary(props_spring, numbNodes = 5L))

## 
## Component sizes
## ```````````````        
## size: 50
##    #:  1
## ______________________________
## Global network properties
## `````````````````````````
## Whole network:
##                                  
## Number of components      1.00000
## Clustering coefficient    0.20010
## Modularity                0.61445
## Positive edge percentage 80.00000
## Edge density              0.06531
## Natural connectivity      0.02556
## Vertex connectivity       1.00000
## Edge connectivity         1.00000
## Average dissimilarity*    0.97871
## Average path length**     3.26313
## -----
## *: Dissimilarity = 1 - edge weight
## **: Path length = Units with average dissimilarity
## 
## ______________________________
## Clusters
## - In the whole network
## - Algorithm: cluster_fast_greedy
## ```````````````````````````````` 
##                   
## name:  1 2  3 4  5
##    #: 14 7 10 8 11
## 
## ______________________________
## Hubs
## - In alphabetical/numerical order
## - Based on empirical quantiles of centralities
## ```````````````````````````````````````````````       
##  193874
##  322860
##  91996 
## 
## ______________________________
## Centrality measures
## - In decreasing order
## - Computed for the complete network
## ````````````````````````````````````
## Degree (unnormalized):
##         
## 91996  9
## 164629 5
## 154360 5
## 218070 5
## 17815  5
## 
## Betweenness centrality (normalized):
##               
## 91996  0.53656
## 25715  0.38350
## 88629  0.33163
## 434721 0.30612
## 434720 0.27891
## 
## Closeness centrality (normalized):
##               
## 91996  0.65420
## 322860 0.55059
## 193874 0.54965
## 158075 0.54771
## 154241 0.53375
## 
## Eigenvector centrality (normalized):
##               
## 91996  1.00000
## 322860 0.68751
## 193874 0.65340
## 154360 0.63810
## 164629 0.62672

pdf(paste0(wdir, "fig/network-spring-16s.pdf"), width = 30, height = 15)
p <- plot(props_spring,
          edgeInvisFilter = "threshold",
          edgeInvisPar = 0.01,
          nodeColor = "cluster",
          nodeSize = "eigenvector",
          repulsion = 0.8,
          rmSingles = TRUE,
          labelScale = FALSE,
          cexLabels = 1.6,
          nodeSizeSpread = 3,
          cexNodes = 2,
          hubBorderCol = "darkgray",
          title1 = "Network on species level with SPRING associations",
          showTitle = TRUE,
          cexTitle = 2.3)
legend(0.57, 1.1, cex = 2.2, title = "estimated association:",
       legend = c("+","-"), lty = 1, lwd = 3, col = c("#009900","red"), 
       bty = "n", horiz = FALSE)
dev.off()

#different transparency value is added to edges with an absolute weight below and above the `cut` value
print(p$q1$Arguments$cut)

biom.16s.genus.top50 <- prune_taxa(names(sort(taxa_sums(biom.16s.genus), decreasing = TRUE)[1:50]), biom.16s.genus)
net_spring_genus <- netConstruct(biom.16s.genus.top50,
                           #filtTax = "highestFreq",
                           #filtTaxPar = list(highestFreq = 50),
                           taxRank = "Genus",
                           measure = "spring",
                           measurePar = list(nlambda=50, 
                                             rep.num=50,
                                             Rmethod = "approx"),
                           normMethod = "none", 
                           zeroMethod = "none",
                           sparsMethod = "none", 
                           dissFunc = "signed",
                           verbose = 3,
                           seed = 123456)

## Checking input arguments ... Done.
## 50 taxa and 24 samples remaining.
## 
## Calculate 'spring' associations ... 
## Done.

#Adjacency values computed using the "signed" transformation (values different from 0 and 1 will be edges in the network):
hist(net_spring_genus$assoMat1, 100, xlim = c(-1, 1), ylim = c(0, 400),
     xlab = "Estimated correlation", 
     main = "Estimated correlations after sparsification")

net_spring_genus_unsigned <- netConstruct(data = net_spring_genus$assoEst1,
                           dataType = "correlation",
                           sparsMethod = "threshold",
                           thresh = 0.3,
                           dissFunc = "unsigned",
                           verbose = 3)

## Checking input arguments ... Done.
## 
## Sparsify associations via 'threshold' ... Done.

#Adjacency values computed using the "unsigned" transformation:
hist(net_spring_genus_unsigned$adjaMat1, 100, ylim = c(0, 400),
     xlab = "Adjacency values", 
     main = "Adjacencies (with \"unsigned\" transformation)")

temp.nw <- net_spring_genus
props_spring_genus <- netAnalyze(temp.nw, 
                           centrLCC = TRUE,
                           clustMethod = "cluster_fast_greedy",
                           #hubPar = "eigenvector",
                           hubPar = "betweenness",
                           hubQuant = 0.9,
                           weightDeg = FALSE, normDeg = FALSE)

#?summary.microNetProps
print(summary(props_spring_genus, numbNodes = 5L))

## 
## Component sizes
## ```````````````        
## size: 50
##    #:  1
## ______________________________
## Global network properties
## `````````````````````````
## Whole network:
##                                  
## Number of components      1.00000
## Clustering coefficient    0.22046
## Modularity                0.61130
## Positive edge percentage 79.72973
## Edge density              0.06041
## Natural connectivity      0.02528
## Vertex connectivity       1.00000
## Edge connectivity         1.00000
## Average dissimilarity*    0.98042
## Average path length**     2.94850
## -----
## *: Dissimilarity = 1 - edge weight
## **: Path length = Units with average dissimilarity
## 
## ______________________________
## Clusters
## - In the whole network
## - Algorithm: cluster_fast_greedy
## ```````````````````````````````` 
##                      
## name:  1  2 3 4 5 6 7
##    #: 12 11 6 4 6 5 6
## 
## ______________________________
## Hubs
## - In alphabetical/numerical order
## - Based on empirical quantiles of centralities
## ```````````````````````````````````````````````                        
##  Leptolyngbya_LEGE-06070
##  Marinovum              
##  Phaeobacter            
##  Pseudarcobacter        
##  RS62_marine_group      
## 
## ______________________________
## Centrality measures
## - In decreasing order
## - Computed for the complete network
## ````````````````````````````````````
## Degree (unnormalized):
##                          
## Marinovum               6
## Leptolyngbya_LEGE-06070 6
## Sulfurovum              5
## Pseudarcobacter         5
## Phaeobacter             5
## 
## Betweenness centrality (normalized):
##                                
## Pseudarcobacter         0.23044
## RS62_marine_group       0.22364
## Leptolyngbya_LEGE-06070 0.21259
## Marinovum               0.21173
## Phaeobacter             0.19728
## 
## Closeness centrality (normalized):
##                                
## Marinovum               0.57124
## Phaeobacter             0.54792
## Pseudarcobacter         0.54395
## RS62_marine_group       0.53124
## Leptolyngbya_LEGE-06070 0.52883
## 
## Eigenvector centrality (normalized):
##                          
## Sulfurovum        1.00000
## RS62_marine_group 0.97937
## NS5_marine_group  0.82058
## Thiomicrorhabdus  0.76638
## Polaribacter      0.72628

pdf(paste0(wdir, "fig/network-spring-genus-16s.pdf"), width = 30, height = 15)
p <- plot(props_spring_genus,
          edgeInvisFilter = "threshold",
          edgeInvisPar = 0.01,
          nodeColor = "cluster",
          colorVec =  c("#804000", "#56B4E9", "#009960", "#999999", '#ff7f00', "#CC79A7", "#E69F00"),
          nodeSize = "eigenvector",
          repulsion = 0.8,
          rmSingles = TRUE,
          labelScale = FALSE,
          cexLabels = 2.4,
          nodeSizeSpread = 3,
          cexNodes = 2,
          hubBorderCol = "darkgray",
          #title1 = "Network on genus level with Pearson correlations",
          showTitle = FALSE,
          #cexTitle = 2.3,
          posCol = "#0072B2",
          negCol = '#e41a1c')
legend(0.75, 1.1, cex = 2.4, title = "estimated association:",
       legend = c("+","-"), lty = 1, lwd = 3, col = c("#0072B2", '#e41a1c'), 
       bty = "n", horiz = FALSE)
dev.off()

## png 
##   2

#different transparency value is added to edges with an absolute weight below and above the `cut` value
print(p$q1$Arguments$cut)

##      75% 
## 0.344691

knitr::kable(as.matrix(sort(colSums(net_spring_genus$normCounts1), decreasing = TRUE)[1:20]))

Marivita	123729.200
Aureispira	80162.459
Lentibacter	51138.237
Geitlerinema_PCC-7105	31473.064
Methylophaga	30384.441
NS3a_marine_group	22247.805
Methylomicrobium	19798.463
Flavobacterium	16975.246
Saprospiraceae_uncultured	16357.548
Lewinella	14769.526
Cryomorphaceae_uncultured	12163.325
Marivivens	9780.107
Planktomarina	9761.932
Thalassospira	8487.833
Lutimaribacter	8360.902
Winogradskyella	8272.679
Cyanobium_PCC-6307	7024.799
Meridianimaribacter	6564.303
Symphothece_PCC-7002	6339.630
Pseudorhodobacter	5181.541

SPRING for time point 11 hours

#https://github.com/stefpeschel/NetCoMi
library(NetCoMi)
biom.16s.genus.top50.11hr <- prune_samples(sample_data(biom.16s.genus.top50) %>% as.data.frame() %>% as.matrix() %>% as.data.frame() %>% filter(Time_point == 11) %>% rownames(), biom.16s.genus.top50)

net_spring_genus.11 <- netConstruct(biom.16s.genus.top50.11hr,
                           #filtTax = "highestFreq",
                           #filtTaxPar = list(highestFreq = 50),
                           taxRank = "Genus",
                           measure = "spring",
                           measurePar = list(nlambda=50, 
                                             rep.num=50,
                                             Rmethod = "approx"),
                           normMethod = "none", 
                           zeroMethod = "none",
                           sparsMethod = "none", 
                           dissFunc = "signed",
                           verbose = 3,
                           seed = 123456)

## Checking input arguments ... Done.
## 50 taxa and 8 samples remaining.
## 
## Calculate 'spring' associations ... 
## Done.

#Adjacency values computed using the "signed" transformation (values different from 0 and 1 will be edges in the network):
hist(net_spring_genus.11$assoMat1, 100, xlim = c(-1, 1), ylim = c(0, 400),
     xlab = "Estimated correlation", 
     main = "Estimated correlations after sparsification")

temp.nw <- net_spring_genus.11
props_spring_genus.11 <- netAnalyze(temp.nw, 
                           centrLCC = TRUE,
                           clustMethod = "cluster_fast_greedy",
                           hubPar = "eigenvector",
                           weightDeg = FALSE, normDeg = FALSE)

#?summary.microNetProps
print(summary(props_spring_genus.11, numbNodes = 5L))

## 
## Component sizes
## ```````````````                   
## size: 17 6 5 4 2  1
##    #:  1 1 1 1 3 12
## ______________________________
## Global network properties
## `````````````````````````
## Largest connected component (LCC):
##                                  
## Relative LCC size         0.34000
## Clustering coefficient    0.30749
## Modularity                0.43878
## Positive edge percentage 90.47619
## Edge density              0.15441
## Natural connectivity      0.07955
## Vertex connectivity       1.00000
## Edge connectivity         1.00000
## Average dissimilarity*    0.95218
## Average path length**     2.48950
## 
## Whole network:
##                                  
## Number of components     19.00000
## Clustering coefficient    0.36825
## Modularity                0.72875
## Positive edge percentage 95.00000
## Edge density              0.03265
## Natural connectivity      0.02355
## -----
## *: Dissimilarity = 1 - edge weight
## **: Path length = Units with average dissimilarity
## 
## ______________________________
## Clusters
## - In the whole network
## - Algorithm: cluster_fast_greedy
## ```````````````````````````````` 
##                           
## name:  0 1 2 3 4 5 6 7 8 9
##    #: 12 7 7 6 5 4 2 3 2 2
## 
## ______________________________
## Hubs
## - In alphabetical/numerical order
## - Based on empirical quantiles of centralities
## ```````````````````````````````````````````````               
##  Phaeobacter   
##  Psychroserpens
##  Tropicibacter 
## 
## ______________________________
## Centrality measures
## - In decreasing order
## - Centrality of disconnected components is zero
## ````````````````````````````````````````````````
## Degree (unnormalized):
##                   
## Phaeobacter      5
## Tropicibacter    5
## Leucothrix       4
## Psychroserpens   4
## Rhodobacteraceae 3
## 
## Betweenness centrality (normalized):
##                                     
## Phaeobacter                  0.64167
## Leucothrix                   0.56667
## Methylophaga                 0.32500
## Planktomarina                0.23333
## Flavobacteriaceae_uncultured 0.23333
## 
## Closeness centrality (normalized):
##                                     
## Phaeobacter                  0.81422
## Tropicibacter                0.81414
## Leucothrix                   0.74242
## Psychroserpens               0.67666
## Flavobacteriaceae_uncultured 0.64401
## 
## Eigenvector centrality (normalized):
##                                     
## Phaeobacter                  1.00000
## Tropicibacter                0.99982
## Psychroserpens               0.67955
## Leucothrix                   0.65893
## Flavobacteriaceae_uncultured 0.59254

pdf(paste0(wdir, "fig/network-spring-genus-11hr-16s.pdf"), width = 30, height = 15)
p <- plot(props_spring_genus.11,
          edgeInvisFilter = "threshold",
          edgeInvisPar = 0.01,
          nodeColor = "cluster",
          colorVec =  c("#804000", "#56B4E9", "#009960", "#999999", "#357090", "#CC79A7", "#E69F00"),
          nodeSize = "eigenvector",
          repulsion = 0.8,
          rmSingles = TRUE,
          labelScale = FALSE,
          cexLabels = 1.6,
          nodeSizeSpread = 3,
          cexNodes = 2,
          hubBorderCol = "darkgray",
          #title1 = "Network on genus level with Pearson correlations",
          showTitle = FALSE,
          #cexTitle = 2.3,
          posCol = "#0072B2",
          negCol = '#e41a1c')
legend(0.75, 1.1, cex = 2.2, title = "estimated association:",
       legend = c("+","-"), lty = 1, lwd = 3, col = c("#0072B2", '#e41a1c'), 
       bty = "n", horiz = FALSE)
dev.off()

## png 
##   2

#different transparency value is added to edges with an absolute weight below and above the `cut` value
print(p$q1$Arguments$cut)

##       75% 
## 0.3308134

knitr::kable(as.matrix(sort(colSums(net_spring_genus.11$normCounts1), decreasing = TRUE)[1:20]))

Marivita	93758.648
Aureispira	25102.082
Methylophaga	17800.597
Meridianimaribacter	5905.413
Winogradskyella	5088.945
Methylomicrobium	4118.112
Leucothrix	3891.838
Lutimaribacter	3275.711
Lewinella	3131.341
Psychroserpens	3079.585
Taeseokella	2880.727
Cellulophaga	2872.356
Flavobacterium	2466.225
Cyanobium_PCC-6307	2459.497
Rhodobacteraceae_uncultured	2316.024
Marinovum	2315.527
Rhodobacteraceae	2002.619
Hoeflea	1963.873
Methylobacter	1924.185
Phaeobacter	1905.210

Pearson’s correlation

#While with the “signed” transformation, positive correlated taxa are likely to belong to the same cluster, with the “unsigned” transformation clusters contain strongly positive and negative correlated taxa.
net_pears <- netConstruct(biom.16s,  
                          filtTax = "highestFreq",
                          filtTaxPar = list(highestFreq = 50),
                          measure = "pearson",
                          normMethod = "clr",
                          zeroMethod = "multRepl",
                          sparsMethod = "threshold",
                          thresh = 0.3,
                          dissFunc = "signed",
                          verbose = 3)

## Checking input arguments ... Done.
## Data filtering ...
## 928 taxa removed.
## 50 taxa and 24 samples remaining.
## 
## Zero treatment:
## Execute multRepl() ... Done.
## 
## Normalization:
## Execute clr(){SpiecEasi} ... Done.
## 
## Calculate 'pearson' associations ... Done.
## 
## Sparsify associations via 'threshold' ... Done.

#Adjacency values computed using the "signed" transformation (values different from 0 and 1 will be edges in the network):
hist(net_pears$assoMat1, 100, xlim = c(-1, 1), ylim = c(0, 400),
     xlab = "Estimated correlation", 
     main = "Estimated correlations after sparsification")

net_pears_unsigned <- netConstruct(data = net_pears$assoEst1,
                           dataType = "correlation",
                           sparsMethod = "threshold",
                           thresh = 0.3,
                           dissFunc = "unsigned",
                           verbose = 3)

## Checking input arguments ... Done.
## 
## Sparsify associations via 'threshold' ... Done.

#Adjacency values computed using the "unsigned" transformation:
hist(net_pears_unsigned$adjaMat1, 100, ylim = c(0, 400),
     xlab = "Adjacency values", 
     main = "Adjacencies (with \"unsigned\" transformation)")

temp.nw <- net_pears
props_pears <- netAnalyze(temp.nw, 
                           centrLCC = TRUE,
                           clustMethod = "cluster_fast_greedy",
                           hubPar = "eigenvector",
                           weightDeg = FALSE, normDeg = FALSE
                          )

#?summary.microNetProps
print(summary(props_pears, numbNodes = 5L))

## 
## Component sizes
## ```````````````        
## size: 50
##    #:  1
## ______________________________
## Global network properties
## `````````````````````````
## Whole network:
##                                  
## Number of components      1.00000
## Clustering coefficient    0.78484
## Modularity                0.04831
## Positive edge percentage 44.33962
## Edge density              0.69224
## Natural connectivity      0.19868
## Vertex connectivity      13.00000
## Edge connectivity        13.00000
## Average dissimilarity*    0.77277
## Average path length**     0.99918
## -----
## *: Dissimilarity = 1 - edge weight
## **: Path length = Units with average dissimilarity
## 
## ______________________________
## Clusters
## - In the whole network
## - Algorithm: cluster_fast_greedy
## ```````````````````````````````` 
##               
## name:  1  2  3
##    #: 14 20 16
## 
## ______________________________
## Hubs
## - In alphabetical/numerical order
## - Based on empirical quantiles of centralities
## ```````````````````````````````````````````````       
##  154281
##  193874
##  91970 
## 
## ______________________________
## Centrality measures
## - In decreasing order
## - Computed for the complete network
## ````````````````````````````````````
## Degree (unnormalized):
##          
## 154281 43
## 193874 43
## 91970  42
## 205285 42
## 154150 41
## 
## Betweenness centrality (normalized):
##               
## 431434 0.03316
## 193874 0.03316
## 391654 0.02211
## 332216 0.02126
## 154281 0.01956
## 
## Closeness centrality (normalized):
##               
## 154360 1.50977
## 193874 1.50318
## 205285 1.49145
## 91970  1.46304
## 154281 1.45782
## 
## Eigenvector centrality (normalized):
##               
## 193874 1.00000
## 154281 0.99297
## 91970  0.98020
## 205285 0.97734
## 154360 0.96529

pdf(paste0(wdir, "fig/network-pearson-16s.pdf"), width = 30, height = 15)
p <- plot(props_pears,
          edgeInvisFilter = "threshold",
          edgeInvisPar = 0.001,
          nodeColor = "cluster",
          nodeSize = "eigenvector",
          repulsion = 0.8,
          rmSingles = TRUE,
          labelScale = FALSE,
          cexLabels = 1.6,
          nodeSizeSpread = 3,
          cexNodes = 2,
          hubBorderCol = "darkgray",
          title1 = "Network on species level with Pearson correlations",
          showTitle = TRUE,
          cexTitle = 2.3)
legend(0.57, 1.1, cex = 2.2, title = "estimated association:",
       legend = c("+","-"), lty = 1, lwd = 3, col = c("#009900","red"), 
       bty = "n", horiz = FALSE)
dev.off()

#different transparency value is added to edges with an absolute weight below and above the `cut` value
print(p$q1$Arguments$cut)

biom.16s.genus.top50 <- prune_taxa(names(sort(taxa_sums(biom.16s.genus), decreasing = TRUE)[1:50]), biom.16s.genus)
net_pears_genus <- netConstruct(biom.16s.genus.top50,
                          taxRank = "Genus",
                          measure = "pearson",
                          zeroMethod = "multRepl",
                          normMethod = "clr",
                          sparsMethod = "threshold",
                          thresh = 0.3,
                          verbose = 3)

## Checking input arguments ... Done.
## 50 taxa and 24 samples remaining.
## 
## Zero treatment:
## Execute multRepl() ... Done.
## 
## Normalization:
## Execute clr(){SpiecEasi} ... Done.
## 
## Calculate 'pearson' associations ... Done.
## 
## Sparsify associations via 'threshold' ... Done.

#Adjacency values computed using the "signed" transformation (values different from 0 and 1 will be edges in the network):
hist(net_pears_genus$assoMat1, 100, xlim = c(-1, 1), ylim = c(0, 400),
     xlab = "Estimated correlation", 
     main = "Estimated correlations after sparsification")

net_pears_genus_unsigned <- netConstruct(data = net_pears_genus$assoEst1,
                           dataType = "correlation",
                           sparsMethod = "threshold",
                           thresh = 0.3,
                           dissFunc = "unsigned",
                           verbose = 3)

## Checking input arguments ... Done.
## 
## Sparsify associations via 'threshold' ... Done.

#Adjacency values computed using the "unsigned" transformation:
hist(net_pears_genus_unsigned$adjaMat1, 100, ylim = c(0, 400),
     xlab = "Adjacency values", 
     main = "Adjacencies (with \"unsigned\" transformation)")

temp.nw <- net_pears_genus
props_pears_genus <- netAnalyze(temp.nw, 
                           centrLCC = TRUE,
                           clustMethod = "cluster_fast_greedy",
                           hubPar = "eigenvector",
                           weightDeg = FALSE, normDeg = FALSE)

#?summary.microNetProps
print(summary(props_pears_genus, numbNodes = 5L))

## 
## Component sizes
## ```````````````        
## size: 50
##    #:  1
## ______________________________
## Global network properties
## `````````````````````````
## Whole network:
##                                  
## Number of components      1.00000
## Clustering coefficient    0.70118
## Modularity                0.06109
## Positive edge percentage 46.16420
## Edge density              0.60653
## Natural connectivity      0.15972
## Vertex connectivity      11.00000
## Edge connectivity        11.00000
## Average dissimilarity*    0.80366
## Average path length**     1.00555
## -----
## *: Dissimilarity = 1 - edge weight
## **: Path length = Units with average dissimilarity
## 
## ______________________________
## Clusters
## - In the whole network
## - Algorithm: cluster_fast_greedy
## ```````````````````````````````` 
##               
## name:  1  2  3
##    #: 17 18 15
## 
## ______________________________
## Hubs
## - In alphabetical/numerical order
## - Based on empirical quantiles of centralities
## ```````````````````````````````````````````````                  
##  Lentibacter      
##  RS62_marine_group
##  Sulfurovum       
## 
## ______________________________
## Centrality measures
## - In decreasing order
## - Computed for the complete network
## ````````````````````````````````````
## Degree (unnormalized):
##                     
## Lentibacter       42
## Pseudarcobacter   38
## RS62_marine_group 38
## Sulfurovum        37
## Tropicibacter     36
## 
## Betweenness centrality (normalized):
##                           
## Marinovum          0.02721
## Lentibacter        0.02551
## Marivivens         0.02296
## Cyanobium_PCC-6307 0.01956
## Lewinella          0.01616
## 
## Closeness centrality (normalized):
##                         
## Sulfurovum       1.42923
## NS5_marine_group 1.40426
## Lentibacter      1.36705
## Methylobacter    1.36133
## Tropicibacter    1.33878
## 
## Eigenvector centrality (normalized):
##                          
## Lentibacter       1.00000
## Sulfurovum        0.97226
## RS62_marine_group 0.92932
## NS5_marine_group  0.91158
## Pseudarcobacter   0.90952

pdf(paste0(wdir, "fig/network-pearson-genus-16s.pdf"), width = 30, height = 15)
p <- plot(props_pears_genus,
          edgeInvisFilter = "threshold",
          edgeInvisPar = 0.5,
          nodeColor = "cluster",
          colorVec =  c("#804000", "#357090", "#009960"),
          nodeSize = "eigenvector",
          repulsion = 0.8,
          rmSingles = TRUE,
          labelScale = FALSE,
          cexLabels = 1.6,
          nodeSizeSpread = 3,
          cexNodes = 2,
          hubBorderCol = "darkgray",
          #title1 = "Network on genus level with Pearson correlations",
          showTitle = FALSE,
          #cexTitle = 2.3)
          posCol = "#0072B2",
          negCol = '#e41a1c')
legend(0.75, 1.1, cex = 2.2, title = "estimated association:",
       legend = c("+","-"), lty = 1, lwd = 3, col = c("#0072B2", '#e41a1c'), 
       bty = "n", horiz = FALSE)
dev.off()

## png 
##   2

#different transparency value is added to edges with an absolute weight below and above the `cut` value
print(p$q1$Arguments$cut)

## [1] 0.5890202

knitr::kable(as.matrix(sort(colSums(net_pears_genus$normCounts1), decreasing = TRUE)[1:20]))

Marivita	75.932138
Lentibacter	46.134143
Flavobacterium	34.910881
Lewinella	28.127237
NS3a_marine_group	27.840424
Methylophaga	27.692529
Aureispira	24.007935
Lutimaribacter	18.709552
Geitlerinema_PCC-7105	12.684909
Cyanobium_PCC-6307	12.209718
Saprospiraceae_uncultured	9.137318
Winogradskyella	8.136879
Methylomicrobium	8.059228
Psychroserpens	6.897760
Pseudorhodobacter	5.296707
Cryomorphaceae_uncultured	4.584473
Marivivens	4.465113
Planktomarina	3.711412
Taeseokella	2.791422
Hydrogenophaga	1.265531

Compare genera that were significant in differential abundance analysis

#lightVSocean: 11
my.genus.1 <- res_sig.aldex.lightVSocean %>% dplyr::select(Genus) %>% as.data.frame()
my.genus.1$lightVSocean <- 1
#lightVSdark: 21
my.genus.2 <- tax_table(biom.16s.m.genus)[glm.eff.lightVSdark[glm.eff.lightVSdark$effect >= 1.5,] %>% row.names(),] %>% as.data.frame() %>% dplyr::select(Genus)
my.genus.2$lightVSdark <- 1
#light.time.trend: 21
my.genus.3 <- dge.5.u %>% dplyr::select(Genus) %>% as.data.frame
my.genus.3$light.time.trend <- 1
#merge
my.genus <- merge(x = my.genus.1, y = my.genus.2, by = "Genus", all = TRUE)
my.genus <- merge(x = my.genus, y = my.genus.3, by = "Genus", all = TRUE)

biom.16s.genus.1 <- biom.16s.genus
sample_data(biom.16s.genus.1) <- sample_data(biom.16s.genus) %>% as.data.frame() %>% as.matrix() %>% as.data.frame() %>% unite(culture_setup_time, c("Culture_setup", "Time_point"))
temp.genus <- merge_samples(biom.16s.genus.1, "culture_setup_time", fun=sum) %>% transform_sample_counts(function(x) {x / sum(x)}) %>% psmelt()
temp.abund <- data.frame()
for(i in 1:length(my.genus$Genus)) {
  temp.out <- temp.genus %>% dplyr::filter(Genus %in% my.genus$Genus[i]) %>% dplyr::select(Sample,Abundance,Genus)
  temp.abund <- rbind(temp.abund, temp.out)
}
my.genus.abund <- reshape2::acast(temp.abund, Genus~Sample, value.var = "Abundance")[,c(11,6,7,4,5,9,10,8,2,3,1)] %>% as.data.frame() %>% tibble::rownames_to_column(var = "Genus")
my.genus.light.16s <- merge(x = my.genus, y = my.genus.abund, by = "Genus", all = TRUE)
knitr::kable(as.matrix(my.genus.light.16s))

Genus	lightVSocean	lightVSdark	light.time.trend	ocean_0	light+clean_gas_5	light+clean_gas_7	light+clean_gas_11	light+clean_gas_14	light+used_gas_5	light+used_gas_7	light+used_gas_11	dark+clean_gas_5	dark+clean_gas_7	dark+clean_gas_11
Ahrensia	NA	NA	1	0.000000e+00	0.0000000000	0.0002755937	0.0012485268	0.000000000	0.0000000000	0.0000000000	0.0017333274	0.0000000000	0.000000000	0.0000000000
Aphanizomenon_NIES81	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0024296440	0.000000000	0.0000000000	0.0000000000	0.0028496620	0.0000000000	0.000000000	0.0000000000
Aureispira	1	1	NA	6.808685e-03	0.0212254783	0.3910692863	0.4739771827	0.000000000	0.0176424624	0.4763975828	0.2968543447	0.0000000000	0.004895320	0.0000000000
Cohaesibacter	NA	NA	1	0.000000e+00	0.0000000000	0.0006749333	0.0030096130	0.000000000	0.0000000000	0.0000000000	0.0023237684	0.0000000000	0.000000000	0.0003355920
Congregibacter	NA	1	NA	0.000000e+00	0.0000000000	0.0002322213	0.0011860258	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000
Cyanobium_PCC-6307	NA	1	NA	2.090557e-02	0.0042945492	0.0049360049	0.0652519415	0.042124248	0.0040058679	0.0026274961	0.0111591612	0.0000000000	0.000000000	0.0001446189
Cyclobacterium	NA	NA	1	8.361058e-05	0.0000000000	0.0008338526	0.0090586584	0.010778457	0.0000000000	0.0000000000	0.0064543753	0.0000000000	0.000000000	0.0018638068
Denitromonas	NA	NA	1	0.000000e+00	0.0000000000	0.0009830893	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0043392431	0.0019635142	0.000000000	0.0015089648
Devosia	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0004476291	0.000000000	0.0000000000	0.0000000000	0.0011458209	0.0000000000	0.000000000	0.0002455133
Flexithrix	NA	1	NA	0.000000e+00	0.0036042295	0.0032408892	0.0000000000	0.000000000	0.0005576115	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000
Geitlerinema_PCC-7105	1	1	NA	2.856894e-03	0.4369261401	0.0466434451	0.0156534819	0.013997676	0.3475910717	0.0300407531	0.0179048376	0.0006570818	0.000000000	0.0000000000
Geminocystis_PCC-6308	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0020425959	0.000000000	0.0000000000	0.0000000000	0.0010466521	0.0000000000	0.000000000	0.0000000000
Gimesia	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0047030673	0.000000000	0.0000000000	0.0000000000	0.0033184462	0.0000000000	0.000000000	0.0000000000
Hoeflea	NA	NA	1	6.799129e-04	0.0015239639	0.0042782162	0.0192128650	0.064235061	0.0014487109	0.0021087878	0.0207822045	0.0018202109	0.002825818	0.0040161491
Labrenzia	NA	NA	1	1.552768e-04	0.0000000000	0.0017825184	0.0030800261	0.000000000	0.0000000000	0.0017992925	0.0040907596	0.0008018876	0.000000000	0.0027660224
Leptolyngbya_LEGE-06070	1	1	NA	2.216495e-04	0.0279353024	0.0049286984	0.0022635154	0.000000000	0.0285922643	0.0040799291	0.0013624441	0.0000000000	0.000000000	0.0000000000
Marivita	NA	NA	1	1.466339e-02	0.0976943417	0.1303832950	0.2146588184	0.299145251	0.1089507411	0.1430977866	0.4013861934	0.0308147473	0.079643158	0.4551946777
Membranicola	NA	NA	1	0.000000e+00	0.0000000000	0.0001763800	0.0019654980	0.000000000	0.0000000000	0.0005459890	0.0034007996	0.0000000000	0.000000000	0.0015698152
Methylomicrobium	NA	NA	1	4.107774e-03	0.0000000000	0.0003282668	0.0012333458	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.1986096679	0.250173474	0.0253570848
Methylophaga	NA	NA	1	3.379498e-03	0.0007062692	0.0050298790	0.0084394570	0.082254218	0.0008083664	0.0012557918	0.0035625787	0.1631236464	0.129891161	0.1082383874
Nautella	NA	NA	1	0.000000e+00	0.0000000000	0.0002162934	0.0000000000	0.000000000	0.0000000000	0.0008171678	0.0013259973	0.0000000000	0.000000000	0.0028750247
Nioella	NA	NA	1	4.563681e-04	0.0007750110	0.0032633822	0.0057504545	0.024132413	0.0017131879	0.0016280435	0.0046061824	0.0037895812	0.003820766	0.0035431194
NS2b_marine_group	NA	1	NA	0.000000e+00	0.0000000000	0.0022541692	0.0003727375	0.021066518	0.0006067598	0.0016204349	0.0003967816	0.0000000000	0.000000000	0.0000000000
Porphyrobacter	1	1	1	0.000000e+00	0.0000000000	0.0001899139	0.0023693512	0.000000000	0.0000000000	0.0012566816	0.0000000000	0.0000000000	0.000000000	0.0000000000
Pseudohoeflea	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0010752881	0.000000000	0.0000000000	0.0000000000	0.0007783249	0.0000000000	0.000000000	0.0001511979
Pseudorhodobacter	1	NA	NA	9.518934e-04	0.0073002743	0.0235046912	0.0081642848	0.021447560	0.0082964378	0.0345157831	0.0068608816	0.0031892845	0.003693128	0.0021981458
Rheinheimera	NA	NA	1	2.114121e-04	0.0000000000	0.0005473819	0.0012469665	0.011358610	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0001135262
Schizothrix_LEGE_07164	1	1	NA	0.000000e+00	0.0007606409	0.0009768852	0.0000000000	0.000000000	0.0011680427	0.0005776241	0.0000000000	0.0000000000	0.000000000	0.0000000000
Symphothece_PCC-7002	1	1	NA	4.971324e-04	0.0558599164	0.0226521751	0.0023310233	0.000000000	0.0576151105	0.0084967991	0.0059078838	0.0000000000	0.000000000	0.0000000000
Synechocystis_PCC-6803	1	1	1	0.000000e+00	0.0000000000	0.0000000000	0.0015656554	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000
Taeseokella	NA	NA	1	4.677293e-04	0.0026763916	0.0047277154	0.0328121404	0.009217898	0.0011426820	0.0070170485	0.0367128905	0.0011930301	0.004288548	0.0036569506
Thalassospira	1	1	NA	7.347337e-04	0.0020072938	0.0471862906	0.0134617628	0.000000000	0.0029291942	0.0678548179	0.0310070252	0.0000000000	0.000000000	0.0000000000

#darkVSocean: 17
my.genus.4 <- res_sig.aldex.darkVSocean %>% dplyr::select(Genus) %>% as.data.frame()
my.genus.4$darkVSocean <- 1
#darkVSlight: 51
my.genus.5 <- tax_table(biom.16s.m.genus)[glm.eff.lightVSdark[glm.eff.lightVSdark$effect <= -1.5,] %>% row.names(),] %>% as.data.frame() %>% dplyr::select(Genus)
my.genus.5$darkVSlight <- 1
#dark.time.trend: 49
my.genus.6 <- dge.6.u %>% dplyr::select(Genus) %>% as.data.frame
my.genus.6$dark.time.trend <- 1
#merge
my.genus <- merge(x = my.genus.4, y = my.genus.5, by = "Genus", all = TRUE)
my.genus <- merge(x = my.genus, y = my.genus.6, by = "Genus", all = TRUE)

temp.abund <- data.frame()
for(i in 1:length(my.genus$Genus)) {
  temp.out <- temp.genus %>% dplyr::filter(Genus %in% my.genus$Genus[i]) %>% dplyr::select(Sample,Abundance,Genus)
  temp.abund <- rbind(temp.abund, temp.out)
}
my.genus.abund <- reshape2::acast(temp.abund, Genus~Sample, value.var = "Abundance")[,c(11,6,7,4,5,9,10,8,2,3,1)] %>% as.data.frame() %>% tibble::rownames_to_column(var = "Genus")
my.genus.dark.16s <- merge(x = my.genus, y = my.genus.abund, by = "Genus", all = TRUE)
knitr::kable(as.matrix(my.genus.dark.16s))

Genus	darkVSocean	darkVSlight	dark.time.trend	ocean_0	light+clean_gas_5	light+clean_gas_7	light+clean_gas_11	light+clean_gas_14	light+used_gas_5	light+used_gas_7	light+used_gas_11	dark+clean_gas_5	dark+clean_gas_7	dark+clean_gas_11
Aequorivita	1	1	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	1.832269e-03
Albimonas	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0003998426	0.0000000000	0.0018285727	0.001345524	2.944809e-03
Alcanivorax	1	NA	NA	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0038442545	0.002124497	2.011676e-03
Antarctobacter	1	1	1	8.116697e-05	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0028121981	0.000000000	2.009017e-03
Aquibacter	NA	1	NA	3.156040e-03	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0019088860	0.000000000	4.895079e-04
Arenibacter	NA	NA	1	1.163180e-03	0.0000000000	0.0012916132	0.0000000000	0.020757124	0.0000000000	0.0000000000	0.0000000000	0.0011389587	0.003186910	1.095420e-02
Bizionia	NA	NA	1	5.079168e-04	0.0000000000	0.0000000000	0.0000000000	0.005515630	0.0000000000	0.0000000000	0.0010516700	0.0005828828	0.004748518	4.758394e-03
Candidatus_Megaira	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0005704324	0.000000000	1.978612e-03
Cellulophaga	NA	NA	1	1.075482e-03	0.0007806143	0.0011316863	0.0004697176	0.000000000	0.0004173870	0.0013700723	0.0014665556	0.0036618264	0.010797998	1.744993e-02
Cohaesibacter	NA	NA	1	0.000000e+00	0.0000000000	0.0006749333	0.0030096130	0.000000000	0.0000000000	0.0000000000	0.0023237684	0.0000000000	0.000000000	3.355920e-04
Cyclobacterium	NA	NA	1	8.361058e-05	0.0000000000	0.0008338526	0.0090586584	0.010778457	0.0000000000	0.0000000000	0.0064543753	0.0000000000	0.000000000	1.863807e-03
Denitromonas	NA	NA	1	0.000000e+00	0.0000000000	0.0009830893	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0043392431	0.0019635142	0.000000000	1.508965e-03
Flavobacteriaceae	1	1	1	1.465568e-03	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0022198774	0.005559746	4.333731e-03
Hoeflea	NA	NA	1	6.799129e-04	0.0015239639	0.0042782162	0.0192128650	0.064235061	0.0014487109	0.0021087878	0.0207822045	0.0018202109	0.002825818	4.016149e-03
Labrenzia	NA	NA	1	1.552768e-04	0.0000000000	0.0017825184	0.0030800261	0.000000000	0.0000000000	0.0017992925	0.0040907596	0.0008018876	0.000000000	2.766022e-03
Legionellaceae_uncultured	NA	NA	1	0.000000e+00	0.0027012178	0.0028538941	0.0000000000	0.000000000	0.0009356290	0.0046546499	0.0000000000	0.0002673780	0.000000000	4.439936e-03
Leucothrix	1	1	1	2.672806e-03	0.0008177871	0.0000000000	0.0000000000	0.000000000	0.0006132225	0.0000000000	0.0000000000	0.0064480522	0.012109591	2.419180e-02
Marinicella	NA	NA	1	2.056224e-04	0.0000000000	0.0003129874	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0008715804	0.003337973	5.296964e-03
Mariniflexile	1	1	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0006591556	0.002593059	3.042980e-03
Marinovum	NA	NA	1	5.102455e-04	0.0012782168	0.0031034346	0.0037769213	0.000000000	0.0008687154	0.0044648891	0.0042497784	0.0115157462	0.007952507	1.274805e-02
Marivita	NA	NA	1	1.466339e-02	0.0976943417	0.1303832950	0.2146588184	0.299145251	0.1089507411	0.1430977866	0.4013861934	0.0308147473	0.079643158	4.551947e-01
Membranicola	NA	NA	1	0.000000e+00	0.0000000000	0.0001763800	0.0019654980	0.000000000	0.0000000000	0.0005459890	0.0034007996	0.0000000000	0.000000000	1.569815e-03
Meridianimaribacter	1	1	1	7.487941e-04	0.0013691579	0.0005527849	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0029245680	0.019121129	3.670825e-02
Mesoflavibacter	1	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.007918675	3.666090e-03
Methylobacter	1	1	NA	5.468674e-04	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0438901766	0.040019563	1.196080e-02
Methylomicrobium	1	1	NA	4.107774e-03	0.0000000000	0.0003282668	0.0012333458	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.1986096679	0.250173474	2.535708e-02
Methylomonadaceae_uncultured	1	1	NA	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0082250754	0.013203332	5.923953e-04
Methylomonas	1	1	NA	0.000000e+00	0.0000000000	0.0017159768	0.0006856936	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0150801648	0.015066388	2.634419e-03
Methylophaga	1	1	NA	3.379498e-03	0.0007062692	0.0050298790	0.0084394570	0.082254218	0.0008083664	0.0012557918	0.0035625787	0.1631236464	0.129891161	1.082384e-01
Methylosarcina	1	1	NA	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0053555983	0.005473982	7.487388e-05
Muricauda	NA	NA	1	0.000000e+00	0.0000000000	0.0014466318	0.0000000000	0.016041084	0.0000000000	0.0005947820	0.0000000000	0.0000000000	0.000000000	1.619021e-03
Nannocystaceae_uncultured	1	1	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0006582078	0.001632063	1.110616e-03
Nannocystis	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	7.748460e-04
Nautella	NA	NA	1	0.000000e+00	0.0000000000	0.0002162934	0.0000000000	0.000000000	0.0000000000	0.0008171678	0.0013259973	0.0000000000	0.000000000	2.875025e-03
Neptunomonas	NA	NA	1	1.353942e-04	0.0000000000	0.0008513155	0.0005442763	0.000000000	0.0017000879	0.0043217679	0.0106348511	0.0006595328	0.001338421	1.644501e-03
Nioella	NA	NA	1	4.563681e-04	0.0007750110	0.0032633822	0.0057504545	0.024132413	0.0017131879	0.0016280435	0.0046061824	0.0037895812	0.003820766	3.543119e-03
Oceanobacterium	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0006021070	0.000000000	1.341588e-03
Oricola	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0005442818	0.000000000	0.0000000000	0.0000000000	0.0012449576	0.0000000000	0.000000000	3.357265e-04
Phaeobacter	NA	NA	1	2.009532e-04	0.0000000000	0.0004764591	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0027059497	0.000000000	1.184285e-02
Primorskyibacter	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0005834410	0.000000000	1.903747e-03
Pseudohoeflea	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0010752881	0.000000000	0.0000000000	0.0000000000	0.0007783249	0.0000000000	0.000000000	1.511979e-04
Psychroserpens	NA	NA	1	1.392088e-03	0.0066090015	0.0055332905	0.0013989326	0.005940084	0.0040755943	0.0059378838	0.0020921003	0.0029410488	0.010420519	1.842288e-02
Rhodobacteraceae	NA	NA	1	5.648481e-04	0.0000000000	0.0022997705	0.0010142864	0.006512690	0.0007569106	0.0044767990	0.0041107715	0.0055404839	0.000000000	1.137303e-02
Rhodobacteraceae_uncultured	NA	NA	1	6.301788e-03	0.0000000000	0.0017093934	0.0015060674	0.000000000	0.0000000000	0.0003141622	0.0003524818	0.0071847670	0.000000000	1.402671e-02
Roseobacter	NA	NA	1	3.498325e-04	0.0000000000	0.0006497742	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0015646925	0.003364529	1.785765e-03
Ruegeria	NA	NA	1	0.000000e+00	0.0000000000	0.0004328169	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0007281610	0.000000000	4.203980e-03
Sagittula	NA	NA	1	1.786031e-04	0.0000000000	0.0008082698	0.0010889112	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0029867921	0.005385985	3.560316e-03
Schleiferia	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	6.434068e-04
Sedimentitalea	1	NA	1	0.000000e+00	0.0000000000	0.0002101179	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0011499838	0.000000000	9.641517e-04
SM1A02	NA	NA	1	0.000000e+00	0.0000000000	0.0006803159	0.0000000000	0.016875810	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	1.204017e-03
Taeseokella	NA	NA	1	4.677293e-04	0.0026763916	0.0047277154	0.0328121404	0.009217898	0.0011426820	0.0070170485	0.0367128905	0.0011930301	0.004288548	3.656951e-03
Tenacibaculum	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0003503334	0.002321306	1.059759e-03
Tropicibacter	1	1	NA	6.143363e-04	0.0000000000	0.0001705410	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0124296965	0.010003630	5.452494e-03
Vitellibacter	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	8.207555e-04
Winogradskyella	NA	NA	1	2.137090e-03	0.0012266796	0.0232761781	0.0005511586	0.012842383	0.0008125185	0.0043722932	0.0013221244	0.0030721625	0.017461028	3.124321e-02
Xanthomarina	NA	NA	1	0.000000e+00	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0000000000	0.000000000	4.259341e-04
Youngimonas	1	1	NA	1.522092e-04	0.0000000000	0.0000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000	0.0000000000	0.0111975727	0.004940597	4.408582e-03

List of all genera with relative abundances per culture setup and time point

biom.16s.genus.1 <- biom.16s.genus
sample_data(biom.16s.genus.1) <- sample_data(biom.16s.genus) %>% as.data.frame() %>% as.matrix() %>% as.data.frame() %>% unite(culture_setup_time, c("Culture_setup", "Time_point"))
temp.genus <- merge_samples(biom.16s.genus.1, "culture_setup_time", fun=sum) %>% transform_sample_counts(function(x) {x / sum(x)}) %>% psmelt()
temp.abund.all <- temp.genus %>% dplyr::select(Sample,Abundance,Genus)
rel.abund.16s <- reshape2::acast(temp.abund.all, Genus~Sample, value.var = "Abundance")
write.table(as.matrix(rel.abund.16s), file = paste0(wdir, "tab-rel-abund-culture-time-16s.csv"), col.names = T, row.names = T, sep = ",", quote = F)

Relative abundance of selected genera

#relative abundance of sequences from methanotrophic bacteria (27--33\%) in the middle of the experiment (day 5--7)
sample_data(biom.16s)$culture_setup_time_point <- mapply(paste0, as.character(get_variable(biom.16s, "Culture_setup")), as.character(get_variable(biom.16s, "Time_point")), collapse = "_")
temp.mob <- merge_samples(biom.16s, "culture_setup_time_point") %>% tax_glom(taxrank = "Phylum") %>% transform_sample_counts(function(x) {x / sum(x)}) %>% psmelt() %>% filter(Phylum %in% "Methanotroph") %>% dplyr::select(Sample, Abundance, Phylum)

#At the end of the dark incubation, other chemoheterotrophic microbes (especially the Proteobacteria \textit{Marivita} and \textit{Methylophaga}) had outcompeted the methanotrophs.
temp.dark.11 <- merge_samples(biom.16s, "culture_setup_time_point") %>% tax_glom(taxrank = "Genus") %>% transform_sample_counts(function(x) {x / sum(x)}) %>% psmelt() %>% filter(Sample %in% "dark+clean_gas11" & Abundance>0.01) %>% dplyr::select(Sample, Abundance, Phylum, Genus)

temp.cyano <- merge_samples(biom.16s, "culture_setup_time_point") %>% tax_glom(taxrank = "Phylum") %>% transform_sample_counts(function(x) {x / sum(x)}) %>% psmelt() %>% filter(Phylum %in% "Cyanobacteria") %>% dplyr::select(Sample, Abundance, Phylum)

temp.light.5 <- merge_samples(biom.16s, "culture_setup_time_point") %>% tax_glom(taxrank = "Genus") %>% transform_sample_counts(function(x) {x / sum(x)}) %>% psmelt() %>% filter(Sample %in% "light+clean_gas5" & Abundance>0.01) %>% dplyr::select(Sample, Abundance, Phylum, Genus)

temp.light.11 <- merge_samples(biom.16s, "culture_setup_time_point") %>% tax_glom(taxrank = "Genus") %>% transform_sample_counts(function(x) {x / sum(x)}) %>% psmelt() %>% filter(Sample %in% "light+clean_gas11" & Abundance>0.01) %>% dplyr::select(Sample, Abundance, Phylum, Genus)

my.genus <- c(
  "Geitlerinema_PCC-7105",
  "Methylomicrobium",
  "Methylobacter",
  "Methylomonas",
  "Marivivens",
  "Winogradskyella",
  "Symphothece_PCC-7002",
  "Cyanobium_PCC-6307",
  "Alteromonas"
)
temp.genus <- merge_samples(biom.16s, "Culture_setup") %>% tax_glom(taxrank = "Genus") %>% transform_sample_counts(function(x) {x / sum(x)}) %>% psmelt()
temp.abund <- data.frame()
for(i in 1:length(my.genus)) {
  temp.out <- temp.genus %>% filter(Genus==my.genus[i]) %>%dplyr::select(Sample, Abundance, Genus)
  temp.abund <- rbind(temp.abund, temp.out)
}
knitr::kable(as.matrix(reshape2::acast(temp.abund, Genus~Sample, value.var = "Abundance")))

	dark+clean_gas	light+clean_gas	light+used_gas	ocean
Alteromonas	0.0000000	0.0008821	0.0067583	0.0001008
Cyanobium_PCC-6307	0.0001001	0.0179715	0.0058653	0.0209056
Geitlerinema_PCC-7105	0.0001616	0.1087390	0.1348442	0.0028569
Marivivens	0.0276675	0.0002971	0.0001893	0.0189411
Methylobacter	0.0215455	0.0000000	0.0000000	0.0005469
Methylomicrobium	0.0818416	0.0004492	0.0000000	0.0041078
Methylomonas	0.0064627	0.0011652	0.0000000	0.0000000
Symphothece_PCC-7002	0.0000000	0.0238959	0.0244769	0.0004971
Winogradskyella	0.0234637	0.0144901	0.0021645	0.0021371

#save workspace
save.image(paste0(devdir, "/phyloseq-16s.Rdata"))

session information

sessionInfo()

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3;  LAPACK version 3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Copenhagen
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] splines   stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] NetCoMi_1.1.0               SpiecEasi_1.1.3            
##  [3] DirichletMultinomial_1.46.0 ALDEx2_1.36.0              
##  [5] latticeExtra_0.6-30         zCompositions_1.5.0-3      
##  [7] truncnorm_1.0-9             NADA_1.6-1.1               
##  [9] survival_3.6-4              MASS_7.3-60.2              
## [11] microbiome_1.26.0           factoextra_1.0.7           
## [13] FactoMineR_2.11             ggbiplot_0.6.2             
## [15] knitr_1.47                  ggpubr_0.6.0               
## [17] patchwork_1.2.0             gplots_3.1.3.1             
## [19] edgeR_4.2.0                 limma_3.60.2               
## [21] DESeq2_1.44.0               SummarizedExperiment_1.34.0
## [23] Biobase_2.64.0              MatrixGenerics_1.16.0      
## [25] matrixStats_1.3.0           GenomicRanges_1.56.0       
## [27] GenomeInfoDb_1.40.1         IRanges_2.38.0             
## [29] S4Vectors_0.42.0            BiocGenerics_0.50.0        
## [31] vegan_2.6-6.1               lattice_0.22-6             
## [33] permute_0.9-7               ggrepel_0.9.5              
## [35] RColorBrewer_1.1-3          psadd_0.1.3                
## [37] reshape2_1.4.4              phyloseq_1.48.0            
## [39] lubridate_1.9.3             forcats_1.0.0              
## [41] stringr_1.5.1               dplyr_1.1.4                
## [43] purrr_1.0.2                 readr_2.1.5                
## [45] tidyr_1.3.1                 tibble_3.2.1               
## [47] tidyverse_2.0.0             ggplot2_3.5.1              
## 
## loaded via a namespace (and not attached):
##   [1] nnet_7.3-19             DT_0.33                 Biostrings_2.72.1      
##   [4] TH.data_1.1-2           vctrs_0.6.5             digest_0.6.35          
##   [7] png_0.1-8               corpcor_1.6.10          shape_1.4.6.1          
##  [10] pcaPP_2.0-4             registry_0.5-1          corrplot_0.92          
##  [13] deldir_2.0-4            parallelly_1.37.1       SPRING_1.0.4           
##  [16] foreach_1.5.2           withr_3.0.0             psych_2.4.3            
##  [19] xfun_0.44               doRNG_1.8.6             memoise_2.0.1          
##  [22] emmeans_1.10.2          latentcor_2.0.1         zoo_1.8-12             
##  [25] gtools_3.9.5            pbapply_1.7-2           Formula_1.2-5          
##  [28] KEGGREST_1.44.0         scatterplot3d_0.3-44    httr_1.4.7             
##  [31] rstatix_0.7.2           globals_0.16.3          rhdf5filters_1.16.0    
##  [34] rhdf5_2.48.0            rstudioapi_0.16.0       UCSC.utils_1.0.0       
##  [37] generics_0.1.3          base64enc_0.1-3         zlibbioc_1.50.0        
##  [40] ca_0.71.1               doSNOW_1.0.20           RcppZiggurat_0.1.6     
##  [43] GenomeInfoDbData_1.2.12 quadprog_1.5-8          SparseArray_1.4.8      
##  [46] xtable_1.8-4            ade4_1.7-22             doParallel_1.0.17      
##  [49] evaluate_0.23           S4Arrays_1.4.1          Rfast_2.1.0            
##  [52] preprocessCore_1.66.0   hms_1.1.3               glmnet_4.1-8           
##  [55] pulsar_0.3.11           irlba_2.3.5.1           colorspace_2.1-0       
##  [58] magrittr_2.0.3          viridis_0.6.5           future.apply_1.11.2    
##  [61] cowplot_1.1.3           Hmisc_5.1-3             pillar_1.9.0           
##  [64] nlme_3.1-164            huge_1.3.5              iterators_1.0.14       
##  [67] caTools_1.18.2          compiler_4.4.1          stringi_1.8.4          
##  [70] biomformat_1.32.0       TSP_1.2-4               dendextend_1.17.1      
##  [73] plyr_1.8.9              crayon_1.5.2            abind_1.4-5            
##  [76] timeSeries_4032.109     sn_2.1.1                locfit_1.5-9.9         
##  [79] bit_4.0.5               rootSolve_1.8.2.4       mixedCCA_1.6.2         
##  [82] sandwich_3.1-0          fastcluster_1.2.6       codetools_0.2-20       
##  [85] multcomp_1.4-25         directlabels_2024.1.21  bslib_0.7.0            
##  [88] plotly_4.10.4           multtest_2.60.0         Rcpp_1.0.12            
##  [91] qgraph_1.9.8            magic_1.6-1             interp_1.1-6           
##  [94] leaps_3.1               blob_1.2.4              utf8_1.2.4             
##  [97] fBasics_4032.96         pbivnorm_0.6.0          listenv_0.9.1          
## [100] checkmate_2.3.1         Rdpack_2.6              ggsignif_0.6.4         
## [103] estimability_1.5.1      lavaan_0.6-17           Matrix_1.7-0           
## [106] statmod_1.5.0           tzdb_0.4.0              pkgconfig_2.0.3        
## [109] tools_4.4.1             cachem_1.1.0            rbibutils_2.2.16       
## [112] RSQLite_2.3.7           viridisLite_0.4.2       DBI_1.2.3              
## [115] numDeriv_2016.8-1.1     impute_1.78.0           fastmap_1.2.0          
## [118] rmarkdown_2.27          scales_1.3.0            grid_4.4.1             
## [121] fMultivar_4031.84       broom_1.0.6             sass_0.4.9             
## [124] coda_0.19-4.1           carData_3.0-5           rpart_4.1.23           
## [127] snow_0.4-4              farver_2.1.2            mgcv_1.9-1             
## [130] yaml_2.3.8              VGAM_1.1-11             spatial_7.3-17         
## [133] foreign_0.8-86          cli_3.6.2               webshot_0.5.5          
## [136] lifecycle_1.0.4         mvtnorm_1.2-5           backports_1.5.0        
## [139] BiocParallel_1.38.0     timechange_0.3.0        gtable_0.3.5           
## [142] parallel_4.4.1          ape_5.8                 jsonlite_1.8.8         
## [145] seriation_1.5.5         bitops_1.0-7            multcompView_0.1-10    
## [148] bit64_4.0.5             assertthat_0.2.1        Rtsne_0.17             
## [151] glasso_1.11             doFuture_1.0.1          heatmaply_1.5.0        
## [154] RcppParallel_5.1.7      jquerylib_0.1.4         highr_0.11             
## [157] timeDate_4032.109       orca_1.1-2              lazyeval_0.2.2         
## [160] dynamicTreeCut_1.63-1   htmltools_0.5.8.1       GO.db_3.19.1           
## [163] glue_1.7.0              XVector_0.44.0          microbenchmark_1.4.10  
## [166] mnormt_2.1.1            jpeg_0.1-10             gridExtra_2.3          
## [169] flashClust_1.01-2       igraph_2.0.3            R6_2.5.1               
## [172] fdrtool_1.2.17          labeling_0.4.3          cluster_2.1.6          
## [175] rngtools_1.5.2          Rhdf5lib_1.26.0         DelayedArray_0.30.1    
## [178] tidyselect_1.2.1        plotrix_3.8-4           WGCNA_1.72-5           
## [181] htmlTable_2.4.2         car_3.1-2               AnnotationDbi_1.66.0   
## [184] future_1.33.2           filematrix_1.3          munsell_0.5.1          
## [187] KernSmooth_2.23-24      geometry_0.4.7          data.table_1.15.4      
## [190] htmlwidgets_1.6.4       rlang_1.1.3             fansi_1.0.6